Matt:LabNotes/2014-7-14

From ZhangLabWiki
Revision as of 03:04, 16 July 2014 by >Mzcai (Created page with "==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==...")
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search

RT Primer mRNA Enrich Sequencing Analysis

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)

Bedtools intersect