Matt:LabNotes/2014-10-31

From ZhangLabWiki
Jump to navigation Jump to search

RT Primer RNA-Seq Analysis[edit]

Align with Tophat2 to hg19[edit]

~/scratch/RTprimer_Analysis$ /home/kunzhang/softwares/tophat-2.0.6.Linux_x86_64/tophat2 -o tophat_hg19unmask_Indx26_RanHex --solexa1.3-quals ~/LTS/Genomes/hg19 s_8_1_Indx26.txt
~/scratch/RTprimer_Analysis$ /home/kunzhang/softwares/tophat-2.0.6.Linux_x86_64/tophat2 -o tophat_hg19unmask_Indx27_dT --solexa1.3-quals ~/LTS/Genomes/hg19 s_8_1_Indx27.txt
~/scratch/RTprimer_Analysis$ /home/kunzhang/softwares/tophat-2.0.6.Linux_x86_64/tophat2 -o tophat_hg19unmask_Indx28_FISSEQRT --solexa1.3-quals ~/LTS/Genomes/hg19 s_8_1_Indx28.txt
~/scratch/RTprimer_Analysis$ /home/kunzhang/softwares/tophat-2.0.6.Linux_x86_64/tophat2 -o tophat_hg19unmask_Indx29_Top48 --solexa1.3-quals ~/LTS/Genomes/hg19 s_8_1_Indx29.txt
  • RanHex
 586553 reads; of these:
 586553 (100.00%) were unpaired; of these:
   251516 (42.88%) aligned 0 times
   196622 (33.52%) aligned exactly 1 time
   138415 (23.60%) aligned >1 times
 57.12% overall alignment rate
  • dT
 432896 reads; of these:
 432896 (100.00%) were unpaired; of these:
   103736 (23.96%) aligned 0 times
   215595 (49.80%) aligned exactly 1 time
   113565 (26.23%) aligned >1 times
 76.04% overall alignment rate
  • FISSEQRT
 529281 reads; of these:
 529281 (100.00%) were unpaired; of these:
   199427 (37.68%) aligned 0 times
   205336 (38.80%) aligned exactly 1 time
   124518 (23.53%) aligned >1 times
 62.32% overall alignment rate
  • Top48
 523107 reads; of these:
 523107 (100.00%) were unpaired; of these:
   165872 (31.71%) aligned 0 times
   229706 (43.91%) aligned exactly 1 time
   127529 (24.38%) aligned >1 times
 68.29% overall alignment rate

Samtools Sort and Index[edit]

 samtools sort tophat_hg19unmask_Indx26_RanHex/accepted_hits.bam mapped_RanHex.sorted
 samtools sort tophat_hg19unmask_Indx27_dT/accepted_hits.bam mapped_dT.sorted
 samtools sort tophat_hg19unmask_Indx28_FISSEQRT/accepted_hits.bam mapped_FISSEQRT.sorted
 samtools sort tophat_hg19unmask_Indx29_Top48/accepted_hits.bam mapped_Top48.sorted
 samtools index mapped_RanHex.sorted.bam mapped_RanHex.sorted.bam.bai
 samtools index mapped_dT.sorted.bam mapped_dT.sorted.bam.bai
 samtools index mapped_FISSEQRT.sorted.bam mapped_FISSEQRT.sorted.bam.bai
 samtools index mapped_Top48.sorted.bam mapped_Top48.sorted.bam.bai
 fetchChromSizes hg19 > ~/Genomes/hg19.chrom.sizes

Visual QC[edit]

 bam2wig.py -i mapped_RanHex.sorted.bam -s /home/mzcai/LTS/Genomes/hg19.chrom.sizes -o mapped_RanHex.sorted
 bam2wig.py -i mapped_dT.sorted.bam -s /home/mzcai/LTS/Genomes/hg19.chrom.sizes -o mapped_dT.sorted
 bam2wig.py -i mapped_FISSEQRT.sorted.bam -s /home/mzcai/LTS/Genomes/hg19.chrom.sizes -o mapped_FISSEQRT.sorted
 bam2wig.py -i mapped_Top48.sorted.bam -s /home/mzcai/LTS/Genomes/hg19.chrom.sizes -o mapped_Top48.sorted


Calculate rRNA Overlap[edit]

hg19_rRNA.bed from UCSC table browser[edit]
 split_bam.py -i mapped_RanHex.sorted.bam -r /home/mzcai/LTS/Genomes/hg19_rRNA.bed -o split_hg19rRNAbed_RanHex
 Total records:                                         398868
 split_hg19rRNAbed_RanHex.in.bam (Reads consumed by input gene list):251593
 split_hg19rRNAbed_RanHex.ex.bam (Reads not consumed by input gene list):147275
 split_hg19rRNAbed_RanHex.junk.bam (qcfailed, unmapped reads):0
 split_bam.py -i mapped_dT.sorted.bam -r /home/mzcai/LTS/Genomes/hg19_rRNA.bed -o split_hg19rRNAbed_dT
 Total records:                                         422589
 split_hg19rRNAbed_dT.in.bam (Reads consumed by input gene list):155754
 split_hg19rRNAbed_dT.ex.bam (Reads not consumed by input gene list):266835
 split_hg19rRNAbed_dT.junk.bam (qcfailed, unmapped reads):0
 split_bam.py -i mapped_FISSEQRT.sorted.bam -r /home/mzcai/LTS/Genomes/hg19_rRNA.bed -o split_hg19rRNAbed_FISSEQRT
 Total records:                                         391847
 split_hg19rRNAbed_FISSEQRT.in.bam (Reads consumed by input gene list):220685
 split_hg19rRNAbed_FISSEQRT.ex.bam (Reads not consumed by input gene list):171162
 split_hg19rRNAbed_FISSEQRT.junk.bam (qcfailed, unmapped reads):0
 split_bam.py -i mapped_Top48.sorted.bam -r /home/mzcai/LTS/Genomes/hg19_rRNA.bed -o split_hg19rRNAbed_Top48
 Total records:                                         422471
 split_hg19rRNAbed_Top48.in.bam (Reads consumed by input gene list):237568
 split_hg19rRNAbed_Top48.ex.bam (Reads not consumed by input gene list):184903
 split_hg19rRNAbed_Top48.junk.bam (qcfailed, unmapped reads):0
Homo_sapiens.GRCh37.75.totalrRNA.chr.bed[edit]
 split_bam.py -i mapped_RanHex.sorted.bam -r /home/mzcai/LTS/Genomes/Homo_sapiens.GRCh37.75.totalrRNA.chr.bed -o split_GRCh37totalrRNAchrbed_RanHex
 Total records:                                         398868
 split_GRCh37totalrRNAchrbed_RanHex.in.bam (Reads consumed by input gene list):47
 split_GRCh37totalrRNAchrbed_RanHex.ex.bam (Reads not consumed by input gene list):398821
 split_GRCh37totalrRNAchrbed_RanHex.junk.bam (qcfailed, unmapped reads):0
 split_bam.py -i mapped_dT.sorted.bam -r /home/mzcai/LTS/Genomes/Homo_sapiens.GRCh37.75.totalrRNA.chr.bed -o split_GRCh37totalrRNAchrbed_dT
 Total records:                                         422589
 split_GRCh37totalrRNAchrbed_dT.in.bam (Reads consumed by input gene list):13
 split_GRCh37totalrRNAchrbed_dT.ex.bam (Reads not consumed by input gene list):422576
 split_GRCh37totalrRNAchrbed_dT.junk.bam (qcfailed, unmapped reads):0
 split_bam.py -i mapped_FISSEQRT.sorted.bam -r /home/mzcai/LTS/Genomes/Homo_sapiens.GRCh37.75.totalrRNA.chr.bed -o split_GRCh37totalrRNAchrbed_FISSEQRT
 Total records:                                         391847
 split_GRCh37totalrRNAchrbed_FISSEQRT.in.bam (Reads consumed by input gene list):43
 split_GRCh37totalrRNAchrbed_FISSEQRT.ex.bam (Reads not consumed by input gene list):391804
 split_GRCh37totalrRNAchrbed_FISSEQRT.junk.bam (qcfailed, unmapped reads):0
 split_bam.py -i mapped_Top48.sorted.bam -r /home/mzcai/LTS/Genomes/Homo_sapiens.GRCh37.75.totalrRNA.chr.bed -o split_GRCh37totalrRNAchrbed_Top48
 Total records:                                         422471
 split_GRCh37totalrRNAchrbed_Top48.in.bam (Reads consumed by input gene list):16
 split_GRCh37totalrRNAchrbed_Top48.ex.bam (Reads not consumed by input gene list):422455
 split_GRCh37totalrRNAchrbed_Top48.junk.bam (qcfailed, unmapped reads):0

Map with Tophat2 to hg19 --report-secondary-alignments[edit]

  • rRNA sequences are often in repeat regions and so will have multiple alignments
    • try reporting all (up to 20) alignments
 /home/kunzhang/softwares/tophat-2.0.6.Linux_x86_64/tophat2 --report-secondary-alignments -o tophat_2ndalign_hg19unmask_Indx26_RanHex --solexa1.3-quals ~/LTS/Genomes/hg19 s_8_1_Indx26.txt
 /home/kunzhang/softwares/tophat-2.0.6.Linux_x86_64/tophat2 --report-secondary-alignments -o tophat_2ndalign_hg19unmask_Indx27_dT --solexa1.3-quals ~/LTS/Genomes/hg19 s_8_1_Indx27.txt
 /home/kunzhang/softwares/tophat-2.0.6.Linux_x86_64/tophat2 --report-secondary-alignments -o tophat_2ndalign_hg19unmask_Indx28_FISSEQRT --solexa1.3-quals ~/LTS/Genomes/hg19 s_8_1_Indx28.txt
 /home/kunzhang/softwares/tophat-2.0.6.Linux_x86_64/tophat2 --report-secondary-alignments -o tophat_2ndalign_hg19unmask_Indx29_Top48 --solexa1.3-quals ~/LTS/Genomes/hg19 s_8_1_Indx29.txt
 samtools sort tophat_2ndalign_hg19unmask_Indx26_RanHex/accepted_hits.bam mapped_2ndalign_RanHex.sorted
 samtools sort tophat_2ndalign_hg19unmask_Indx27_dT/accepted_hits.bam mapped_2ndalign_dT.sorted
 samtools sort tophat_2ndalign_hg19unmask_Indx28_FISSEQRT/accepted_hits.bam mapped_2ndalign_FISSEQRT.sorted
 samtools sort tophat_2ndalign_hg19unmask_Indx29_Top48/accepted_hits.bam mapped_2ndalign_Top48.sorted
 samtools index mapped_2ndalign_RanHex.sorted.bam mapped_2ndalign_RanHex.sorted.bam.bai
 samtools index mapped_2ndalign_dT.sorted.bam mapped_2ndalign_dT.sorted.bam.bai
 samtools index mapped_2ndalign_FISSEQRT.sorted.bam mapped_2ndalign_FISSEQRT.sorted.bam.bai
 samtools index mapped_2ndalign_Top48.sorted.bam mapped_2ndalign_Top48.sorted.bam.bai

Calculate rRNA Overlap[edit]

hg19_rRNA.bed from UCSC table browser[edit]
 split_bam.py -i mapped_2ndalign_RanHex.sorted.bam -r /home/mzcai/LTS/Genomes/hg19_rRNA.bed -o split_hg19rRNAbed_2ndalign_RanHex
 Total records:                                         617797
 split_hg19rRNAbed_2ndalign_RanHex.in.bam (Reads consumed by input gene list):372244
 split_hg19rRNAbed_2ndalign_RanHex.ex.bam (Reads not consumed by input gene list):245553
 split_hg19rRNAbed_2ndalign_RanHex.junk.bam (qcfailed, unmapped reads):0
 split_bam.py -i mapped_2ndalign_dT.sorted.bam -r /home/mzcai/LTS/Genomes/hg19_rRNA.bed -o split_hg19rRNAbed_2ndalign_dT
 Total records:                                         674492
 split_hg19rRNAbed_2ndalign_dT.in.bam (Reads consumed by input gene list):225058
 split_hg19rRNAbed_2ndalign_dT.ex.bam (Reads not consumed by input gene list):449434
 split_hg19rRNAbed_2ndalign_dT.junk.bam (qcfailed, unmapped reads):0
 split_bam.py -i mapped_2ndalign_FISSEQRT.sorted.bam -r /home/mzcai/LTS/Genomes/hg19_rRNA.bed -o split_hg19rRNAbed_2ndalign_FISSEQRT
 Total records:                                         578621
 split_hg19rRNAbed_2ndalign_FISSEQRT.in.bam (Reads consumed by input gene list):312791
 split_hg19rRNAbed_2ndalign_FISSEQRT.ex.bam (Reads not consumed by input gene list):265830
 split_hg19rRNAbed_2ndalign_FISSEQRT.junk.bam (qcfailed, unmapped reads):0
 split_bam.py -i mapped_2ndalign_Top48.sorted.bam -r /home/mzcai/LTS/Genomes/hg19_rRNA.bed -o split_hg19rRNAbed_2ndalign_Top48
 Total records:                                         633213
 split_hg19rRNAbed_2ndalign_Top48.in.bam (Reads consumed by input gene list):351066
 split_hg19rRNAbed_2ndalign_Top48.ex.bam (Reads not consumed by input gene list):282147
 split_hg19rRNAbed_2ndalign_Top48.junk.bam (qcfailed, unmapped reads):0
Homo_sapiens.GRCh37.75.totalrRNA.chr.bed[edit]
 split_bam.py -i mapped_2ndalign_RanHex.sorted.bam -r /home/mzcai/LTS/Genomes/Homo_sapiens.GRCh37.75.totalrRNA.chr.bed -o split_GRCh37totalrRNAchrbed_2ndalign_RanHex
 Total records:                                         617797
 split_GRCh37totalrRNAchrbed_2ndalign_RanHex.in.bam (Reads consumed by input gene list):47
 split_GRCh37totalrRNAchrbed_2ndalign_RanHex.ex.bam (Reads not consumed by input gene list):617750
 split_GRCh37totalrRNAchrbed_2ndalign_RanHex.junk.bam (qcfailed, unmapped reads):0

 split_bam.py -i mapped_2ndalign_dT.sorted.bam -r /home/mzcai/LTS/Genomes/Homo_sapiens.GRCh37.75.totalrRNA.chr.bed -o split_GRCh37totalrRNAchrbed_2ndalign_dT
 Total records:                                         674492
 split_GRCh37totalrRNAchrbed_2ndalign_dT.in.bam (Reads consumed by input gene list):13
 split_GRCh37totalrRNAchrbed_2ndalign_dT.ex.bam (Reads not consumed by input gene list):674479
 split_GRCh37totalrRNAchrbed_2ndalign_dT.junk.bam (qcfailed, unmapped reads):0
 split_bam.py -i mapped_2ndalign_FISSEQRT.sorted.bam -r /home/mzcai/LTS/Genomes/Homo_sapiens.GRCh37.75.totalrRNA.chr.bed -o split_GRCh37totalrRNAchrbed_2ndalign_FISSEQRT
 Total records:                                         578621
 split_GRCh37totalrRNAchrbed_2ndalign_FISSEQRT.in.bam (Reads consumed by input gene list):43
 split_GRCh37totalrRNAchrbed_2ndalign_FISSEQRT.ex.bam (Reads not consumed by input gene list):578578
 split_GRCh37totalrRNAchrbed_2ndalign_FISSEQRT.junk.bam (qcfailed, unmapped reads):0
 split_bam.py -i mapped_2ndalign_Top48.sorted.bam -r /home/mzcai/LTS/Genomes/Homo_sapiens.GRCh37.75.totalrRNA.chr.bed -o split_GRCh37totalrRNAchrbed_2ndalign_Top48
 Total records:                                         633213
 split_GRCh37totalrRNAchrbed_2ndalign_Top48.in.bam (Reads consumed by input gene list):17
 split_GRCh37totalrRNAchrbed_2ndalign_Top48.ex.bam (Reads not consumed by input gene list):633196
 split_GRCh37totalrRNAchrbed_2ndalign_Top48.junk.bam (qcfailed, unmapped reads):0

Overlap with intersectBed[edit]

  • Compare with other results
 ~/softwares/bedtools-2.20.1/bin/intersectBed -abam mapped_2ndalign_RanHex.sorted.bam -b ~/Genomes/Homo_sapiens.GRCh37.75.totalrRNA.chr.gtf > intersect_tophat2ndalign_rRNAchr_RanHex.bed -bed
 ~/softwares/bedtools-2.20.1/bin/intersectBed -abam mapped_2ndalign_dT.sorted.bam -b ~/Genomes/Homo_sapiens.GRCh37.75.totalrRNA.chr.gtf > intersect_tophat2ndalign_rRNAchr_dT.bed -bed
 ~/softwares/bedtools-2.20.1/bin/intersectBed -abam mapped_2ndalign_FISSEQRT.sorted.bam -b ~/Genomes/Homo_sapiens.GRCh37.75.totalrRNA.chr.gtf > intersect_tophat2ndalign_rRNAchr_FISSEQRT.bed -bed  
 ~/softwares/bedtools-2.20.1/bin/intersectBed -abam mapped_2ndalign_Top48.sorted.bam -b ~/Genomes/Homo_sapiens.GRCh37.75.totalrRNA.chr.gtf > intersect_tophat2ndalign_rRNAchr_Top48.bed -bed

Results[edit]

File:RTprimerAnalysis RSeQCResults.PNG

Check method by aligning to hg19.masked[edit]

  • Very few rRNA reads should map to hg19.masked because repeat regions (usually containing rRNA genes) are masked
 /home/kunzhang/softwares/tophat-2.0.6.Linux_x86_64/tophat2 -o tophat_hg19mask_Indx26_RanHex --solexa1.3-quals ~/LTS/Genomes/hg19.masked s_8_1_Indx26.txt
 samtools sort tophat_hg19mask_Indx26_RanHex/accepted_hits.bam mapped_maskedRanHex.sorted
 samtools index mapped_maskedRanHex.sorted.bam mapped_maskedRanHex.sorted.bam.bai
 split_bam.py -i mapped_maskedRanHex.sorted.bam -r /home/mzcai/LTS/Genomes/hg19_rRNA.bed -o split_hg19rRNAbed_maskedRanHex
 Total records:                                         86356
 split_hg19rRNAbed_maskedRanHex.in.bam (Reads consumed by input gene list):0
 split_hg19rRNAbed_maskedRanHex.ex.bam (Reads not consumed by input gene list):86356
 split_hg19rRNAbed_maskedRanHex.junk.bam (qcfailed, unmapped reads):0
  • As expected there were 0 rRNA
    • Validates hg19_rRNA.bed