Data Analysis Hands-on

docker run -i -t -v /path/to/tutorial:/tutorial -w /tutorial ceciliaklein/teaching:uvic

(main) root@516a2307ed2b:/tutorial#

export PATH=$PATH:/tutorial/teaching-utils/;

more /tutorial/refs/gencode.vM4.gtf

grep Sec23a /tutorial/refs/gencode.vM4.gtf > Sec23a.gtf

awk '$3=="transcript"' Sec23a.gtf | wc -l

awk '$3=="exon"' Sec23a.gtf | wc -l

mkdir -p fastqc

fastqc -o fastqc -f fastq /tutorial/data/mouse_brain_rep1_1.fastq.gz

fastqc/mouse_brain_rep1_1_fastqc.html

mkdir alignments

mkdir genomeIndex

STAR --runThreadN 2 \
     --runMode genomeGenerate \
     --genomeDir genomeIndex \
     --genomeFastaFiles /tutorial/refs/genome.fa \
     --sjdbGTFfile /tutorial/refs/anno.gtf \
     --genomeSAindexNbases 11

STAR --runThreadN 2 \
     --genomeDir genomeIndex \
     --readFilesIn /tutorial/data/mouse_brain_rep1_1.fastq.gz \
                   /tutorial/data/mouse_brain_rep1_2.fastq.gz \
     --outSAMunmapped Within \
     --outFilterType BySJout \
     --outSAMattributes NH HI AS NM MD \
     --readFilesCommand pigz -p2 -dc \
     --outSAMtype BAM SortedByCoordinate \
     --quantMode TranscriptomeSAM \
     --outFileNamePrefix alignments/mouse_brain_rep1_

samtools view -h alignments/mouse_brain_rep1_Aligned.sortedByCoord.out.bam | more

cat alignments/mouse_brain_rep1_Log.final.out

                                 Started job on |	Dec 27 16:01:03
                             Started mapping on |	Dec 27 16:01:04
                                    Finished on |	Dec 27 16:02:16
       Mapping speed, Million of reads per hour |	5.00

                          Number of input reads |	100000
                      Average input read length |	202
                                    UNIQUE READS:
                   Uniquely mapped reads number |	90162
                        Uniquely mapped reads % |	90.16%
                          Average mapped length |	200.46
                       Number of splices: Total |	37153
            Number of splices: Annotated (sjdb) |	1851
                       Number of splices: GT/AG |	36794
                       Number of splices: GC/AG |	190
                       Number of splices: AT/AC |	2
               Number of splices: Non-canonical |	167
                      Mismatch rate per base, % |	0.23%
                         Deletion rate per base |	0.01%
                        Deletion average length |	1.59
                        Insertion rate per base |	0.00%
                       Insertion average length |	1.39
                             MULTI-MAPPING READS:
        Number of reads mapped to multiple loci |	1713
             % of reads mapped to multiple loci |	1.71%
        Number of reads mapped to too many loci |	100
             % of reads mapped to too many loci |	0.10%
                                  UNMAPPED READS:
       % of reads unmapped: too many mismatches |	0.00%
                 % of reads unmapped: too short |	7.96%
                     % of reads unmapped: other |	0.06%
                                  CHIMERIC READS:
                       Number of chimeric reads |	0
                            % of chimeric reads |	0.00%

STAR --runThreadN 2 \
     --runMode inputAlignmentsFromBAM \
     --inputBAMfile alignments/mouse_brain_rep1_Aligned.sortedByCoord.out.bam \
     --outWigType bedGraph \
     --outWigStrand Unstranded \
     --outFileNamePrefix alignments/mouse_brain_rep1_ \
     --outWigReferencesPrefix chr

bedGraphToBigWig alignments/mouse_brain_rep1_Signal.Unique.str1.out.bg \
                 /tutorial/refs/genome.fa.fai \
                 alignments/mouse_brain_rep1_uniq.bw

genome: Mouse    assembly: Dec 2011 (GRCm38/mm10)    [mm10]

track autoScale=on name=RNAseq,brain,PRNAembryoBrain1 color=204,204,0 viewLimits=0:80 maxHeightPixels=30 type=bigWig visibility=2 bigDataUrl=https://public-docs.crg.es/rguigo/Data/cklein/courses/UVIC/UCSC/PRNAembryoBrain1.bw
track autoScale=on name=RNAseq,brain,PRNAembryoBrain2 color=204,204,0 viewLimits=0:80 maxHeightPixels=30 type=bigWig visibility=2 bigDataUrl=https://public-docs.crg.es/rguigo/Data/cklein/courses/UVIC/UCSC/PRNAembryoBrain2.bw
track autoScale=on name=RNAseq,heart,PRNAembryoHeart1 color=128,0,128 viewLimits=0:80 maxHeightPixels=30 type=bigWig visibility=2 bigDataUrl=https://public-docs.crg.es/rguigo/Data/cklein/courses/UVIC/UCSC/PRNAembryoHeart1.bw
track autoScale=on name=RNAseq,heart,PRNAembryoHeart2 color=128,0,128 viewLimits=0:80 maxHeightPixels=30 type=bigWig visibility=2 bigDataUrl=https://public-docs.crg.es/rguigo/Data/cklein/courses/UVIC/UCSC/PRNAembryoHeart2.bw
track autoScale=on name=RNAseq,liver,PRNAembryoLiver1 color=107,142,35 viewLimits=0:80 maxHeightPixels=30 type=bigWig visibility=2 bigDataUrl=https://public-docs.crg.es/rguigo/Data/cklein/courses/UVIC/UCSC/PRNAembryoLiver1.bw
track autoScale=on name=RNAseq,liver,PRNAembryoLiver2 color=107,142,35 viewLimits=0:80 maxHeightPixels=30 type=bigWig visibility=2 bigDataUrl=https://public-docs.crg.es/rguigo/Data/cklein/courses/UVIC/UCSC/PRNAembryoLiver2.bw

mkdir quantifications

mkdir txIndex

rsem-prepare-reference --gtf /tutorial/refs/anno.gtf \
                      /tutorial/refs/genome.fa \
                      txIndex/RSEMref

rsem-calculate-expression --num-threads 2 \
                          --bam \
                          --paired-end \
                          --estimate-rspd \
                          --forward-prob 0 \
                          --no-bam-output \
                          --seed 12345 \
                          alignments/mouse_brain_rep1_Aligned.toTranscriptome.out.bam \
                          txIndex/RSEMref \
                          quantifications/mouse_brain_rep1

cat /tutorial/data/gene.quantifications.index.txt \
| retrieve_element_rpkms.py --output encode \
                            --organism mouse \
                            --element gene \
                            --value TPM \
                            --output_dir quantifications

cat /tutorial/data/gene.quantifications.index.txt \
| retrieve_element_rpkms.py --output encode \
                            --organism mouse \
                            --element gene \
                            --value expected_count \
                            --output_dir quantifications

# protein-coding
(head -1 ../quantifications/encode.mouse.gene.TPM.idr_NA.tsv; grep -Ff ../protein_coding_IDs.txt ../quantifications/encode.mouse.gene.TPM.idr_NA.tsv) > ../quantifications/encode.mouse.pc.gene.TPM.idr_NA.tsv

# long non-coding
(head -1 ../quantifications/encode.mouse.gene.TPM.idr_NA.tsv; grep -Ff ../lincRNA_IDs.txt ../quantifications/encode.mouse.gene.TPM.idr_NA.tsv) > ../quantifications/encode.mouse.linc.gene.TPM.idr_NA.tsv

rpkm_distribution.R --input_matrix ../quantifications/encode.mouse.gene.TPM.idr_NA.tsv \
                    --log \
                    --pseudocount 0 \
                    --metadata ../data/gene.quantifications.index.tsv \
                    --fill_by tissue \
                    --palette /tutorial/palettes/palTissue.txt \
                    --output all_genes

boxplot.log_T.psd_0.all_genes.pdf

selectMatrixRows.sh /tutorial/protein_coding_IDs.txt /tutorial/quantifications/encode.mouse.gene.TPM.idr_NA.tsv > protein-coding.TPM.tsv

matrix_to_dist.R --input_matrix ../quantifications/encode.mouse.gene.TPM.idr_NA.tsv \
                 --log10 \
                 --cor pearson \
                 --output stdout \
| ggheatmap.R --input_matrix stdin \
              --row_metadata /tutorial/data/gene.quantifications.index.tsv \
              --col_dendro \
              --row_dendro \
              --base_size 10 \
              --matrix_palette /tutorial/palettes/palSequential.txt \
              --rowSide_palette /tutorial/palettes/palTissue.txt \
              --rowSide_by tissue \
              --matrix_fill_limits 0.85,1 \
              --output samples_clustering.pdf

edgeR.analysis.R --input_matrix ../quantifications/encode.mouse.gene.expected_count.idr_NA.tsv \
                 --metadata /tutorial/data/gene.quantifications.index.tsv \
                 --fields tissue \
                 --coefficient 2 \
                 --output brain_X_heart

awk '$NF<0.01 && $2<-10{print $1"\tover_brain_X_heart"}' edgeR.cpm1.n2.brain_X_heart.tsv > edgeR.0.01.over_brain_X_heart.txt

awk '$NF<0.01 && $2>10 {print $1"\tover_heart_X_brain"}' edgeR.cpm1.n2.brain_X_heart.tsv > edgeR.0.01.over_heart_X_brain.txt

wc -l edgeR.0.01.over*.txt

awk '$3=="gene"{ match($0, /gene_id "([^"]+).+gene_type "([^"]+)/, var); print var[1],var[2] }' OFS="\t" /tutorial/refs/gencode.vM4.gtf \
| join.py --file1 stdin \
          --file2 <(cat edgeR.0.01.over*.txt) \
| sed '1igene\tedgeR\tgene_type' > gene.edgeR.tsv

cut -f1 gene.edgeR.tsv \
| tail -n+2 \
| selectMatrixRows.sh - ../quantifications/encode.mouse.gene.TPM.idr_NA.tsv \
| ggheatmap.R --width 5 \
              --height 8 \
              --col_metadata /tutorial/data/gene.quantifications.index.tsv \
              --colSide_by tissue \
              --col_labels labExpId \
              --row_metadata gene.edgeR.tsv \
              --merge_row_mdata_on gene \
              --rowSide_by edgeR,gene_type \
              --row_labels none \
              --log \
              --pseudocount 0.1 \
              --col_dendro \
              --row_dendro \
              --matrix_palette /tutorial/palettes/palDiverging.txt \
              --colSide_palette /tutorial/palettes/palTissue.txt \
              --output heatmap.brain_X_heart.pdf

awk '{split($10,a,/\"|\./); print a[2]}' /tutorial/refs/gencode.vM4.gtf | sort -u > universe.txt

# BP
awk '{split($1,a,"."); print a[1]}' edgeR.0.01.over_brain_X_heart.txt \
| GO_enrichment.R --universe universe.txt \
                  --genes stdin \
                  --categ BP \
                  --output edgeR.over_brain_X_heart \
                  --species mouse

# MF
awk '{split($1,a,"."); print a[1]}' edgeR.0.01.over_brain_X_heart.txt \
| GO_enrichment.R --universe universe.txt \
                  --genes stdin \
                  --categ MF \
                  --output edgeR.over_brain_X_heart \
                  --species mouse

# CC
awk '{split($1,a,"."); print a[1]}' edgeR.0.01.over_brain_X_heart.txt \
| GO_enrichment.R --universe universe.txt \
                  --genes stdin \
                  --categ CC \
                  --output edgeR.over_brain_X_heart \
                  --species mouse

awk 'NR==1{$1="% "$1}{print $1,$2}' edgeR.over_brain_X_heart.BP.tsv

# List of transcript ids
awk '$3=="transcript" && $0~/gene_type "protein_coding"/{ match($0, /transcript_id "([^"]+)/, id); print id[1] }' /tutorial/refs/gencode.vM4.gtf |sort -u > protein_coding_transcript_IDs.txt

# Genome annotation restricted to exon features and filtered by transcript type
cat /tutorial/refs/gencode.vM4.gtf |awk '$3=="exon"' |grep -Ff protein_coding_transcript_IDs.txt > exon-annot.gtf

# Filter transcript TPM matrix
selectMatrixRows.sh protein_coding_transcript_IDs.txt /tutorial/quantifications/encode.mouse.transcript.TPM.idr_NA.tsv > pc-tx.tsv

# Individual transcript expression matrices
for tissue in Brain Heart Liver; do
    selectMatrixColumns.sh PRNAembryo${tissue}1:PRNAembryo${tissue}2 pc-tx.tsv > expr.${tissue}.tsv
done

suppa.py generateEvents -i exon-annot.gtf -e SE SS MX RI FL -o localEvents -f ioe

wc -l localEvents*ioe

event=SE; suppa.py psiPerEvent --total-filter 10 --ioe-file localEvents_${event}_strict.ioe --expression-file pc-tx.tsv -o PSI-${event}

event=SE; for tissue in Brain Heart Liver ;do selectMatrixColumns.sh PRNAembryo${tissue}1:PRNAembryo${tissue}2 PSI-${event}.psi > ${tissue}.${event}.psi;done

event=SE; suppa.py diffSplice --method empirical --input localEvents_${event}_strict.ioe --psi Brain.${event}.psi Heart.${event}.psi Liver.${event}.psi --tpm expr.Brain.tsv expr.Heart.tsv expr.Liver.tsv -c -gc -o DS.${event}

event=SE; awk 'BEGIN{FS=OFS="\t"}NR==1{print }NR>1 && $2!="nan" && ($2>0.7 || $2<-0.7) && $3<0.05{print}' DS.${event}.dpsi
event=SE; awk 'BEGIN{FS=OFS="\t"}NR==1{print }NR>1 && $4!="nan" && ($4>0.7 || $4<-0.7) && $5<0.05{print}' DS.${event}.dpsi
event=SE; awk 'BEGIN{FS=OFS="\t"}NR==1{print }NR>1 && $6!="nan" && ($6>0.7 || $6<-0.7) && $7<0.05{print}' DS.${event}.dpsi

mkdir chip-analysis

cd chip-analysis