Matt:LabNotes/2014-10-31: Difference between revisions
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== | ||
* | *[[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]
- Bed file from gene annotations of Hg19 from Ensembl
- Converted gtf to bed
- Added 'chr' to chromosome names
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]
- Bed file from gene annotations of Hg19 from Ensembl
- Kept only "rRNA" and converted to bed format
- 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[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