Los archivos que ocuparemos en este módulo, los podrás encontrar en un carpeta de drive para su descarga. En ella se encuentran 8 archivos, todos de libre acceso.
Si quieres mayor descripción de los archivos puedes acceder a su fuente mediante
Ver presentación FastQ
Crea una carpeta en la misma raíz en donde está tu carpeta tool por nombre datos.taller copia los archivos fastq que vamos a utilizar.
Muévete a tu carpeta datos.taller:
cd ~/datos.taller
Despliega las primeraz líneas del archivo fastq.
head head_1mill_input_forward_R1_001.fastq
FastQC es el programa más utilizado para revisar la calidad de las secuencias de los experimentos de secuenciación masiva.
Liga: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
Para descargarlo puedes ir a la página de la herramienta FastQC y hacer click en Download Now y escoges la descarga adecuada para tu sistema operativo Linux/Windows o Mac.
FastQC ya está descargado en Drona pero si quieres instalarlo en tu computadora, además de descargarlo, descomprimirlo y tener java instalado en el sistema. Las instrucciones se encuentran en la liga de descarga. Ahí además están las instrucciones de uso.
Hay dos modos de ejecutar FastQC: interactivo y no-interactivo.
El interactivo sólo se puede utilizar si lo tienen instalado en su propia computadora. Por lo que aquí vamos a usarlo de manera no interactiva a través del servidor. La forma no-interactiva es la más común de utilizar porque de esta manera se incluye en los pipelines.
Interactivo:
fastqc
No-interactivo:
fastqc read1 read2
Vamos a revisar la sección de la documentación donde se describen las gráficas (liga) y luego los ejemplos de datos de mala y buena calidad.
Crea una carpeta para guardar los resultados de FastQC:
mkdir salida_fastqc
Ejecuta Fastqc en los archivos que acabas de copiar:
fastqc head_1mill_input_forward_R1_001.fastq head_1mill_input_reverse_R2_001.fastq -o salida_fastqc
Desde tu computadora, puedes revisar los reportes html creados por FastQC en la carpeta salida_fastqc usando un navegador.
Trimmomatic es un programa para remover adaptadores y bases con baja calidad en nuestras secuencias.
Github: https://github.com/usadellab/Trimmomatic
Trimmomatic ya está en drona, pero si quieres descargarlo en tu computadora estos son los pasos:
wget https://github.com/usadellab/Trimmomatic/files/5854859/Trimmomatic-0.39.zip
unzip Trimmomatic-0.39.zip
Copia todo el texto de abajo y guárdalo en un archivo trimm.sh, esto lo puedes hacer con el bloc de notas, NotePad, gedit, TextEdit, nano, o cualquier otro editor de texto plano.
trimmomatic=/home/tools/Trimmomatic-0.39/
read=~/datos.taller/head_1mill_input_forward_R1_001.fastq
mate=~/datos.taller/head_1mill_input_reverse_R2_001.fastq
outdir=~/datos.taller/salida_trimmomatic
mkdir -p ${outdir}
java -jar ${trimmomatic}/trimmomatic-0.39.jar PE \
$read \
$mate \
${outdir}/read_paired.fq \
${outdir}/read_unpaired.fq \
${outdir}/mate_paired.fq \
${outdir}/mate_unpaired.fq \
ILLUMINACLIP:${trimmomatic}/adapters/NexteraPE-PE.fa:2:30:10 \
LEADING:30 \
TRAILING:30 \
MINLEN:75 \
SLIDINGWINDOW:5:20
Si tu archivo lo creaste en tu computadora local, copialo al servidor con el siguiente comando:
Ejecutemos el script que acabamos de crear
cd ~/datos.taller
bash trimm.sh
Revisemos con FastQC como se ven ahora las secuencias
mkdir salida_fastqc_trimmed
fastqc salida_trimmomatic/read_paired.fq salida_trimmomatic/mate_paired.fq -o salida_fastqc_trimmed
De nuevo en tu computadora, revisa los resultados en un navegador y compara con los resultados antes de Trimmomatic.
Es un programa para la cuantificación de datos de bulk o single-cell RNA-seq
Liga: https://pachterlab.github.io/kallisto/
Primero revisa la sección Getting started
Para poder cuantificar Kallisto necesita un índice del transcriptoma de referencia. Kallisto nos proporciona unos ya generados o podemos generar el nuestro
Liga índices: https://github.com/pachterlab/kallisto-transcriptome-indices/releases
Información de los datos que vamos a usar:
Liga: https://rnabio.org/module-01-inputs/0001/05/01/RNAseq_Data/
Ahora si probemos kallisto en una muestra. Para este ejercicio vamos a utlilizar datos de la liga de arriba que ya están descargados.
cd ~/bioinformatics_express
idx="~/bioinformatics_express/chr22.idx"
reads="~/bioinformatics_express/HBR_Rep1_ERCC-Mix2_Build37-ErccTranscripts-chr22.read1.fastq.gz"
mates="~/bioinformatics_express/HBR_Rep1_ERCC-Mix2_Build37-ErccTranscripts-chr22.read2.fastq.gz"
kallisto quant -i ${idx} -o salida_kallisto --rf-stranded $reads $mates
Inspeccionamos lo que salió
ls salida_kallisto
El archivo que nos interesa en este momento es abundance.tsv que lo podemos ahora ver con los comandos more o less.
more salida_kallisto/abundance.tsv
Para hacer lo mismo ahora varias muestras podemos poner el comando de un bucle o loop. Copia el siguiente código a un archivo de texto que se llame hbr_kallisto.sh.
idx="/home/datos.taller/chr22.idx"
for id in HBR_Rep1 HBR_Rep2 HBR_Rep3
do
reads="/home/datos.taller/${id}_ERCC-Mix2_Build37-ErccTranscripts-chr22.read1.fastq.gz"
mates="/home/datos.taller/${id}_ERCC-Mix2_Build37-ErccTranscripts-chr22.read2.fastq.gz"
kallisto quant -i ${idx} -o ~/datos.taller/salida_kallisto_${id} --rf-stranded $reads $mates
done
Ya que tenenos el archivo, desde drona, en datos.sync ejecutamos el siguiente comando:
bash hbr_kallisto.sh
Esto nos creo 3 carpetas con la cuanticación de las muestras que empiezan con HBR.
Ejercicio: ahora intenta haces un archivo uhr_kallisto.sh que haga lo mismo en las muestras que empiezan con UHR.
Descarga la anotación de genes de la página de Ensembl para humanos, achivo gtf, utilizando los comandos que ya conoces obten un archivo que contenga únicamente la columna con el id del transcript, el id del gen, y el nombre del gen para todos los genes en el chromosoma 22. Es importante no tener líneas repetidas.
El archivo una vez descargado es el siguiente Homo_sapiens.GRCh38.111.chr.gtf.gz, dado que termina en gz entonces quiere decir que está comprimido. Para poder verlo podemos usar el comando zcat, que es la variante del comando cat que nos permite ver archivos comprimidos y combinarlo con el comando head para ver las primeras líneas.
zcat Homo_sapiens.GRCh38.111.chr.gtf.gz | head
Para extraer solo las líneas que pertenecen al cromosoma 22 podemos usar grep, y usando ^ garantizamos que las líneas que extraemos son las que empiezan con '22'.
zcat Homo_sapiens.GRCh38.111.chr.gtf.gz | grep '^22'
Ahora, en la tercera columna de nuestro archivo vemos que la información se encuentra a varios niveles por ejemplo CDS, exon, five_prime_utr, gene, start_codon, stop_codon, three_prime_utr y transcript. A nosotros solo nos interesa el nivel de transcript, entonces vamos a extaer solo los renglones que digan la palabra transcript.
zcat Homo_sapiens.GRCh38.111.chr.gtf.gz | grep '^22' | grep -w transcript
Finalmente vemos que la información que queremos esta en la columna 9, y necesitamos dividir esa información en las 3 columnas que necesitamos: transcript_id, gene_id y gene_name.
zcat Homo_sapiens.GRCh38.111.chr.gtf.gz | grep '^22' | grep -w transcript | cut -f9
Ahí obtuvimos la columna 9, pero ahora para dividir la columna 9 en nuevas columnas vemos que la información está separada por ';', entonces podemos hacerlo de nuevo con cut, pero tenemos que decirle a cut que ya no vamos a usar tabs como separador de columnas y vamos a usar ;
zcat Homo_sapiens.GRCh38.111.chr.gtf.gz | grep '^22' | grep -w transcript | cut -f9 | cut -d';' -f1,3,5
Por ahí se colaron algunas líneas que no tienen nombre de gen, entonces para eliminarlas, hacemos lo siguiente
zcat Homo_sapiens.GRCh38.111.chr.gtf.gz | grep '^22' | grep -w transcript | cut -f9 | cut -d';' -f1,3,5 | grep -w 'gene_name'
Finalmente le damos el formato que queremos, de nuevo usando cut
zcat Homo_sapiens.GRCh38.111.chr.gtf.gz | grep '^22' | grep -w transcript | cut -f9 | cut -d';' -f1,3,5 | grep -w 'gene_name' | cut -d' ' -f2,4,6
Ahora borramos las comillas y los ; y reemplazamos los espacios por tabs:
zcat Homo_sapiens.GRCh38.111.chr.gtf.gz | grep '^22' | grep -w transcript | cut -f9 | cut -d';' -f1,3,5 | grep -w 'gene_name' | cut -d' ' -f2,4,6 | sed 's/"//g' | sed 's/;//' | sed 's/\s/\t/g'
Ahora si, ya tenemos lo que queríamos y lo podemos guardar a un archivo
zcat Homo_sapiens.GRCh38.111.chr.gtf.gz | grep '^22' | grep -w transcript | cut -f9 | cut -d';' -f1,3,5 | grep -w 'gene_name' | cut -d' ' -f2,4,6 | sed 's/"//g' | sed 's/;//' | sed 's/\s/\t/g' | sort -u > geneId_transcriptId_geneName_chr22.txt