Matt:LabNotes/2014-10-31: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Mzcai
mNo edit summary
>Mzcai
mNo edit summary
 
Line 1: Line 1:
==RT Primer RNA-Seq Analysis==
==RT Primer RNA-Seq Analysis==
*Library [[Matt:LabNotes/2014-9-25]]
*[[Matt:LabNotes/2014-9-25 | Sequencing Library]]


===Align with Tophat2 to hg19===
===Align with Tophat2 to hg19===
Line 37: Line 37:
   68.29% overall alignment rate
   68.29% overall alignment rate


===Samtools Sort and Index===
====Samtools Sort and Index====
   samtools sort tophat_hg19unmask_Indx26_RanHex/accepted_hits.bam mapped_RanHex.sorted
   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_Indx27_dT/accepted_hits.bam mapped_dT.sorted
Line 50: Line 50:
   fetchChromSizes hg19 > ~/Genomes/hg19.chrom.sizes
   fetchChromSizes hg19 > ~/Genomes/hg19.chrom.sizes


===Visual QC===
====Visual QC====
   bam2wig.py -i mapped_RanHex.sorted.bam -s /home/mzcai/LTS/Genomes/hg19.chrom.sizes -o mapped_RanHex.sorted
   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_dT.sorted.bam -s /home/mzcai/LTS/Genomes/hg19.chrom.sizes -o mapped_dT.sorted
Line 58: Line 58:
<!-- wigToBigWig wigVarStepExample.gz hg19.chrom.sizes myBigWig.bw -->
<!-- wigToBigWig wigVarStepExample.gz hg19.chrom.sizes myBigWig.bw -->


===Calculate rRNA Overlap===
====Calculate rRNA Overlap====


====hg19_rRNA.bed from UCSC table browser====
=====hg19_rRNA.bed from UCSC table browser=====
   split_bam.py -i mapped_RanHex.sorted.bam -r /home/mzcai/LTS/Genomes/hg19_rRNA.bed -o split_hg19rRNAbed_RanHex
   split_bam.py -i mapped_RanHex.sorted.bam -r /home/mzcai/LTS/Genomes/hg19_rRNA.bed -o split_hg19rRNAbed_RanHex
   Total records:                                        398868
   Total records:                                        398868
Line 85: Line 85:
   split_hg19rRNAbed_Top48.junk.bam (qcfailed, unmapped reads):0
   split_hg19rRNAbed_Top48.junk.bam (qcfailed, unmapped reads):0


====Homo_sapiens.GRCh37.75.totalrRNA.chr.bed====
=====Homo_sapiens.GRCh37.75.totalrRNA.chr.bed=====
*Bed file from [[Matt:LabNotes/2014-7-14#Bedtools_intersect | gene annotations of Hg19 from Ensembl]]
*Bed file from [[Matt:LabNotes/2014-7-14#Bedtools_intersect | gene annotations of Hg19 from Ensembl]]
**Converted gtf to bed
**Converted gtf to bed
Line 113: Line 113:
   split_GRCh37totalrRNAchrbed_Top48.junk.bam (qcfailed, unmapped reads):0
   split_GRCh37totalrRNAchrbed_Top48.junk.bam (qcfailed, unmapped reads):0


===Map with Tophat2===
===Map with Tophat2 to hg19 --report-secondary-alignments===
*--report-secondary-alignments
*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====
 
=====hg19_rRNA.bed from UCSC table browser=====
  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=====
*Bed file from [[Matt:LabNotes/2014-7-14#Bedtools_intersect | gene annotations of Hg19 from Ensembl]]
**Kept only "rRNA" and converted to bed format
**[[Matt:LabNotes/2014-7-24#Bedtools_intersect_troubleshooting | Added 'chr' to chromosome names]]
 
  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====
*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===
[[File:RTprimerAnalysis_RSeQCResults.PNG]]
 
===Check method by aligning to hg19.masked===
*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

Latest revision as of 10:06, 3 November 2014

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