Ensamblaje del genoma viral - quipupe/from-assembly-to-nextstrain GitHub Wiki

Descargar el genoma SRR11508492, que fue el primer secuenciado por el Instituto Nacional de Salud del Perú SRR11508492.

prefetch SRR11508492
fasterq-dump.2.10.5 --skip-technical ../SRA/SRR11508492/SRR11508492.sra -e 8 -m 200 -p -3

Ambos programas se encuentran dentro de SRA Toolkit https://github.com/ncbi/sra-tools

Enlace del archivo con los programas compilados para el sistema Ubuntu sratoolkit.

El archivo se descarga como SRR11508492.sra

fasterq-dump separa los archivos en .1 y .2

Modifica el nombre del archivo con el comando mv y luego comprime ambos archivos usando gzip

mv SRR11508492.sra_1.fastq SRR11508492_1.fastq 
mv SRR11508492.sra_2.fastq SRR11508492_2.fastq 

gzip SRR11508492_1.fastq
gzip SRR11508492_2.fastq

Utiliza FastQC para revisar la calidad de las secuencias

fastqc SRR11508492_1.fastq.gz
fastqc SRR11508492_2.fastq.gz

Luego, usa Trimmomatic para hacer el trimming o cortar las secuencias según parámetros de calidad

PE: Paired-end reads

-threads: Numero de cores

LEADING y TRAILING cortan las secuencias en los extremos hasta que encuentren secuencias con calidades mayores a 20

SLIDINGWINDOW utiliza una ventana de 4 bases para calcular si dentro de ellas la calidad promedio es 20, si es menor se corta la secuencia

MINLEN es el tamaño mínimo de las secuencias

java -jar -Xms8g ~/Trimmomatic-0.39/trimmomatic-0.39.jar PE -threads 7 SRR11508492_1.fastq.gz SRR11508492_2.fastq.gz out_forward_PE.fq.gz out_forward_unPE.fq.gz out_reverse_PE.fq.gz out_reverse_unPE.fq.gz ILLUMINACLIP:NexteraPE-PE.fa:2:30:10:2:keepBothReads LEADING:20 TRAILING:20 SLIDINGWINDOW:4:20 MINLEN:30

Obtuve estos resultados

TrimmomaticPE: Started with arguments:
-threads 7 SRR11508492_1.fastq.gz SRR11508492_2.fastq.gz out_forward_PE.fq.gz out_forward_unPE.fq.gz out_reverse_PE.fq.gz out_reverse_unPE.fq.gzILLUMINACLIP:NexteraPE-PE.fa:2:30:10:2:keepBothReads LEADING:20 TRAILING:20 SLIDINGWINDOW:4:20 MINLEN:30
Using PrefixPair: 'AGATGTGTATAAGAGACAG' and 'AGATGTGTATAAGAGACAG'
Using Long Clipping Sequence: 'GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG'
Using Long Clipping Sequence: 'TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG'Using Long Clipping Sequence: 'CTGTCTCTTATACACATCTCCGAGCCCACGAGAC'Using Long Clipping Sequence: 'CTGTCTCTTATACACATCTGACGCTGCCGACGA'
ILLUMINACLIP: Using 1 prefix pairs, 4 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequences																			Quality encoding detected as
phred33
Input Read Pairs: 2359909 Both Surviving: 1578896 (66.90%) Forward Only Surviving: 703350 (29.80%) Reverse Only Surviving: 21051 (0.89%) Dropped: 56612 (2.40%)
TrimmomaticPE: Completed successfully

Vuelve a hacer fastqc en los archivos trimmeados

fastqc out_forward_PE.fq.gz

Bowtie2

Luego, usa bowtie2 para crear un index a partir de la secuencia fasta. El index se llamara virus

bowtie2-build -f file.fasta dbname #ejemplo

bowtie2-build -f wuhan_NC_045512.fa virus

Use como información esta pagina y esta otra para utilizar estos comandos

Luego se procede con el mapeo: -x busca el índice -1 o -2 se refiere a los fastq pareados -U se refiere a los fastq no pareados -S es el archivo de salida en formato sam --no-unal no considerará los reads no alineados -p es el número de procesadores

bowtie2 -x ../reference/virus -1 out_forward_PE.fq.gz -2 out_reverse_PE.fq.gz -U out_forward_unPE.fq.gz,out_reverse_unPE.fq.gz -S ./virus.sam --no-unal -p 7

Obtuve estos resultados

2303297 reads; of these:
1578896 (68.55%) were paired; of these:
1569199 (99.39%) aligned concordantly 0 times
(0.61%) aligned concordantly exactly 1 time
3 (0.00%) aligned concordantly >1 times 
1569199 pairs aligned concordantly 0 times; of these:
(0.02%) aligned discordantly 1 time
1568878 pairs aligned 0 times concordantly or discordantly; of these:
3137756 mates make up the pairs; of these:
(99.99%) aligned 0 times
188 (0.01%) aligned exactly 1 time
140 (0.00%) aligned >1 times
724401 (31.45%) were unpaired; of these:
720372 (99.44%) aligned 0 times
4015 (0.55%) aligned exactly 1 time
14 (0.00%) aligned >1 times
0.63% overall alignment rate

samtools

Usa samtools para

Tranformar el archivo sam a bam
Ordenar el archivo bam
Crear un indice a partir de la secuencia de referencia
Indexar el archivo bam

mv aln-pe.sam virus.sam
samtools view -bS virus.sam > virus.bam #convertir de sam a bam, bam es un formato más ligero 
samtools sort virus.bam -o virus.sorted.bam #el archivo bam es ordenado
samtools faidx virus.fa #crea un índice a partir del genoma de referencia
samtools index virus.sorted.bam #indexa el archivo ordenado

Con samtools también se puede crear un archivo fastq con la información de los reads pareados paired.fq y los reads no pareados single.fq que sí mapearon al genoma de referencia.

samtools fastq -0 /dev/null -s single.fq -N virus.bam > paired.fq

Si necesitas crear archivos separados por cada paired-end entonces usa este comando, no generará los no pareados

samtools/bin/samtools fastq -1 paired1.fq -2 paired2.fq -0 /dev/null -s /dev/null -n virus.sorted.bam
[M::bam2fq_mainloop] discarded 16454 singletons
[M::bam2fq_mainloop] processed 23144 reads

MEGAHIT

También se puede usar MEGAHIT para ensamblar el genoma a partir de los reads mapeados

python3.6 ~/MEGAHIT-1.2.9-Linux-x86_64-static/bin/megahit --12 ./SRR11508492/mejorado/paired.fq -r ./SRR11508492/mejorado/single.fq -o ./assembly2 -t 8

a mí me funcionó este comando, recordar dónde están sus carpetas
python3.6 ../../../MEGAHIT-1.2.9-Linux-x86_64-static/bin/megahit --12 paired.fq -r single.fq -o ./assembly_megahit -t 8

Otros programas, no utilizados

BWA

http://bio-bwa.sourceforge.net/bwa.shtml

#Crea el índice a partir del genoma de referencia
bwa index reference.fasta
#Se crearán estos archivos, que son la database del BWA
reference.fasta.amb  reference.fasta.ann  reference.fasta.bwt  reference.fasta.pac  reference.fasta.sa
#Luego se hace el mapeo, ojo que solo usará los PAIRED reads y no los unpaired, -t controla los procesadores, el proceso tomó 5 minutos el archivo de salida será un bam 
#Este es un intento inicial, pero no lo usaré porque crea un archivo con los reads mapeados como con los no mapeados así que es demasiado pesado
bwa mem -t 6 reference.fasta ../mejorado/out_forward_PE.fq.gz ../mejorado/out_reverse_PE.fq.gz | samtools view -Sb -F 4  > virus.bam

#Este es el comando a usar, -t controla los procesadores. Le estamos dando un parámetro mínimo de calidad y que borre los no mapeados o mal mapeados y que se quede con los mapeados PAIRED.
bwa mem -t 8 reference.fasta ../mejorado/out_forward_PE.fq.gz ../mejorado/out_reverse_PE.fq.gz | samtools view -q 10 -F 1284 -f 0x02 -b > virus.bam 
#Ordenar el archivo de salida (sort)
samtools sort virus.bam -o virus.sorted.bam 
#Los reads mapeados pueden ser exportados como FASTA para ser usados por otro programa
samtools view virus2.sorted.bam | cut -f1,10 | sed 's/^/>/' | sed 's/\t/\n/' > virus.fasta

Para saber más sobre los FLAGS de BWA chequear https://broadinstitute.github.io/picard/explain-flags.html

SPAdes - Ensamblaje de-novo

Usa SPAdes para ensamblar el genoma Necesitarás los archivos paired y unpaired

python3.6 spades.py --12 ../../SRR11508492/reference/paired.fq -s ../../SRR11508492/reference/single.fq -o ../../assembly -t 8

#Version mejorada
python3.6 ../SPAdes-3.14.1/bin/spades.py --isolare --12 ./SRR11508492/mejorado/paired.fq -s ./SRR11508492/mejorado/single.fq -o ./assembly -t 8

QUAST

Para evaluar el assembly se usa este programa, se necesita la secuencia de referencia y el contig/scaffold.

python3.6 /home/quipu/quast-5.0.2/quast.py scaffolds.fasta -r ../reference.fasta

Comparando el ensamblaje de SPAdes

con el de MEGAHIT