Características

Es útil conocer algunas características básicas del ensamble de un genoma, entre ellas el número de contigs, el N50, el tamaño del genoma y la cobertura. Para calcular esto necesitamos el archivo de ensamble (fasta) y uno de los archivos con los que se realizó en ensamble (fastq), para el caso de archivos pair-end.

En todos los comandos a continuación, se podrán usar los archivos del ensamble (aquí llamado FASTAFILE) y uno de los fastq (FASTQFILE).

Número de contigs

$ grep -c ">" FASTAFILE

N50

$ awk '/^>/ {if (seqlen){print seqlen}; print ;seqlen=0;next; } { seqlen += length($0)}END{print seqlen}' FASTAFILE | grep -v ">" | sort -n | awk '{len[i++]=$1;sum+=$1} END {for (j=0;j<i+1;j++) {csum+=len[j]; if (csum>sum/2) {print len[j];break}}}' | sed -r ':L;s=\b([0-9]+)([0-9]{3})\b=\1,\2=g;t L'

Este comando complejo calcula el N50 y lo presenta con comas.

Tamaño del genoma

$ grep -v ">" FASTAFILE | wc | sed 's/ /\t/g' | rev | cut -f1 | rev

Número de secuencias fastq

$ cat FASTQFILE | echo $((`wc -l`/2))

Este comando calcula el número de bases en ambos archivos fastq (.R1.fastq) por duplicado; en caso de querer solo el de uno de los archivos pair-end, modificarlo a:

$ cat FASTQFILE | echo $((`wc -l`/4))

Número de bases en el archivo fastq

$ cat FASTQFILE | paste - - - - | cut -f 2 | tr -d '\n' | wc -c | awk '{print $1*2}'

Este comando calcula el número de bases en el archivo fastq y lo multiplica por dos, estimando que ambos archivos pair-end contienen el mismo número aproximado de bases.

Cobertura

La cobertura de un ensamble se calcula dividiendo el número de bases en ambos archivos fastq entre el número de bases en el ensamble. Usaremos algunos de los cálculos anteriores pasándolos primero a una variable para hacer el cálculo mas eficientemente.

$ TOTAL_FASTQ_BASES=$(cat FASTQFILE | paste - - - - | cut -f 2 | tr -d '\n' | wc -c | awk '{print $1*2}') 

$ TOTAL_CONTIG_BASES=$(grep -v ">" FASTAFILE | wc | sed 's/ /\t/g' | rev | cut -f1 | rev)

$ COVERAGE=$(bc <<< "scale=1; $TOTAL_FASTQ_BASES / $TOTAL_CONTIG_BASES")

$ echo $COVERAGE