Matt:LabNotes/2014-7-14
Jump to navigation
Jump to search
RT Primer mRNA Enrich Sequencing Analysis
- Sequenced libraries of cDNA generated from UHRR with varying RT primers Matt:LabNotes/2014-6-23
Alignment to Reference
Align with Bowtie2 to total rRNA
- Combine rRNA sequences and Mt_rRNA fasta sequences (from biomart, same ones used to design Top48 RT primers for mRNA enrich)
- Build bowtie2 index
mzcai@genome-miner:~/scratch/RanHex_EnrichmRNA$cat mart_export_exons_Mt_rRNA.fa mart_export_exons_rRNA.fa > mart_export_exons_totalrRNA.fa mzcai@genome-miner:~/scratch/RanHex_EnrichmRNA$bowtie2-build mart_export_exons_totalrRNA.fa tot_rRNA
- Bowtie2 alignment
mzcai@genome-miner:~/scratch/RanHex_EnrichmRNA/inVitroRTSeq$ bowtie2 --phred64 -x ~/scratch/RanHex_EnrichmRNA/tot_rRNA -q s_1_1_Indx26.txt > RanHex_bowtie2_totRNA.sam 5441 reads; of these: 5441 (100.00%) were unpaired; of these: 2964 (54.48%) aligned 0 times 2477 (45.52%) aligned exactly 1 time 0 (0.00%) aligned >1 times 45.52% overall alignment rate mzcai@genome-miner:~/scratch/RanHex_EnrichmRNA/inVitroRTSeq$ bowtie2 --phred64 -x ~/scratch/RanHex_EnrichmRNA/tot_rRNA -q s_1_1_Indx27.txt > dT_bowtie2_totRNA.sam 10905 reads; of these: 10905 (100.00%) were unpaired; of these: 9121 (83.64%) aligned 0 times 1783 (16.35%) aligned exactly 1 time 1 (0.01%) aligned >1 times 16.36% overall alignment rate mzcai@genome-miner:~/scratch/RanHex_EnrichmRNA/inVitroRTSeq$ bowtie2 --phred64 -x ~/scratch/RanHex_EnrichmRNA/tot_rRNA -q s_1_1_Indx28.txt > FISSEQRT_bowtie2_totRNA.sam 7024 reads; of these: 7024 (100.00%) were unpaired; of these: 4574 (65.12%) aligned 0 times 2447 (34.84%) aligned exactly 1 time 3 (0.04%) aligned >1 times 34.88% overall alignment rate mzcai@genome-miner:~/scratch/RanHex_EnrichmRNA/inVitroRTSeq$ bowtie2 --phred64 -x ~/scratch/RanHex_EnrichmRNA/tot_rRNA -q s_1_1_Indx29.txt > Top48_bowtie2_totRNA.sam 12158 reads; of these: 12158 (100.00%) were unpaired; of these: 7388 (60.77%) aligned 0 times 4764 (39.18%) aligned exactly 1 time 6 (0.05%) aligned >1 times 39.23% overall alignment rate mzcai@genome-miner:~/scratch/RanHex_EnrichmRNA/inVitroRTSeq$ bowtie2 --phred64 -x ~/scratch/RanHex_EnrichmRNA/tot_rRNA -q s_1_1_Indx32.txt > gDNA_bowtie2_totRNA.sam 111950 reads; of these: 111950 (100.00%) were unpaired; of these: 111132 (99.27%) aligned 0 times 817 (0.73%) aligned exactly 1 time 1 (0.00%) aligned >1 times 0.73% overall alignment rate
- RanHex: 45.52% (5,441 reads)
- dT: 16.36% (10,905 reads)
- FISSEQRT: 34.88% (7,024 reads)
- Top48: 39.23% (12,158 reads)
- gDNA: 0.73% (111,950 reads)
Align with Tophat2 to Hg19
- Build bowtie2 index for Hg19
bowtie2-build ~/TwoGenesPilotExperiment/myHg19/hg19.fa hg19 mv hg19* ~/Genomes
- Tophat 2.0.6 alignment
/home/kunzhang/softwares/tophat-2.0.6.Linux_x86_64/tophat2 --solexa1.3-quals ~/LTS/Genomes/hg19 s_1_1_Indx26.txt /home/kunzhang/softwares/tophat-2.0.6.Linux_x86_64/tophat2 -o tophat_out_Indx27_dT --solexa1.3-quals ~/LTS/Genomes/hg19 s_1_1_Indx27.txt /home/kunzhang/softwares/tophat-2.0.6.Linux_x86_64/tophat2 -o tophat_out_Indx28_FISSEQRT --solexa1.3-quals ~/LTS/Genomes/hg19 s_1_1_Indx28.txt /home/kunzhang/softwares/tophat-2.0.6.Linux_x86_64/tophat2 -o tophat_out_Indx29_Top48 --solexa1.3-quals ~/LTS/Genomes/hg19 s_1_1_Indx29.txt /home/kunzhang/softwares/tophat-2.0.6.Linux_x86_64/tophat2 -o tophat_out_Indx32_gDNA --solexa1.3-quals ~/LTS/Genomes/hg19 s_1_1_Indx32.txt
- Check tophat_out_Indx*_*/logs/bowtie.left_kept_reads.log for alignment rate
- RanHex: 38.9%
- dT: 55.0%
- FISSEQRT: 44.3%
- Top48: 34.0%
- gDNA: 40.2%
Compare Alignments with Genome Annotation
Bioconductor with TranscriptDb from UCSC (hg19, ensGene)
- Use Bioconductor to annotate Tophat accepted_hits.bam with gene name and then count the number of gene names correspond to rRNA
- R script used:
inputfile <- "../tophat_out_Indx32_gDNA/accepted_hits.bam" outputfile <- "../tophat_out_Indx32_gDNA/accepted_ensGene_rawCounts.txt"
library(GenomicFeatures) library(Rsamtools) txdb=makeTranscriptDbFromUCSC(genome='hg19',tablename='ensGene') tx_by_gene=transcriptsBy(txdb,'gene') reads_accepted=readBamGappedAlignments(inputfile) counts_accepted=countOverlaps(tx_by_gene,reads_accepted,ignore.strand=FALSE) toc=data.frame(cDNA=counts_accepted,stringsAsFactors=FALSE) rownames(toc)=names(tx_by_gene) write.table(toc,file=outputfile,sep="\t",col.names=TRUE,row.names=TRUE)