Matt:LabNotes/2014-7-14: Difference between revisions
Jump to navigation
Jump to search
>Mzcai mNo edit summary |
>Mzcai |
||
Line 86: | Line 86: | ||
write.table(toc,file=outputfile,sep="\t",col.names=TRUE,row.names=TRUE) | write.table(toc,file=outputfile,sep="\t",col.names=TRUE,row.names=TRUE) | ||
*Used bash to edit files a little and then excel vlookup to count number of rRNA genes | *Used bash to edit files a little and then excel vlookup to count number of rRNA genes | ||
**[[ | **[[Media:Seq_Tophat_rRNA_Analysis.xlsx | Excel workbook here]] | ||
***First sheet is counts of each gene for ignore.strand=FALSE | ***First sheet is counts of each gene for ignore.strand=FALSE | ||
***Second sheet is counts of each for ignore.strand=TRUE | ***Second sheet is counts of each for ignore.strand=TRUE |
Revision as of 03:37, 16 July 2014
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:
- Tried both ignore.strand=TRUE and ignore.strand=FALSE
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)
- Used bash to edit files a little and then excel vlookup to count number of rRNA genes
- Excel workbook here
- First sheet is counts of each gene for ignore.strand=FALSE
- Second sheet is counts of each for ignore.strand=TRUE
- Third sheet is counts of each rRNA gene for every sequenced sample
- Excel workbook here
Results Summary
- 0 rRNA genes found for every sample
- This result is unlikely to be correct
- TranscriptDb form UCSC missing 26 of the 370 rRNA genes found when I used biomart to get Ensembl genes to design the RT primers
- It's missing many more genes compared to biomart when we look at all genes
- Total number of annotated alignments:
' | Total aligned | Annotated with ignore.strand=FALSE | Annotated with ignore.strand=TRUE |
RanHex | 2116 | 1221 | 2599 |
dT | 5997 | 3964 | 8044 |
FISSEQRT | 3111 | 1913 | 3745 |
Top48 | 4133 | 2555 | 4975 |
gDNA | 44982 | 15441 | 31096 |
Bedtools intersect
- For gene annotations of Hg19 downloaded GTF file from Ensembl FTP:
- Homo_sapiens.GRCh37.75.gtf.gz
- Saved in ~/Genomes/ and gunzipped to ~/Genomes/Homo_sapiens.GRCh37.75.gtf
- Made subset of rRNA only gene annotations
grep rRNA ~/Genomes/Homo_sapiens.GRCh37.75.gtf > ~/Genomes/Homo_sapiens.GRCh37.75.totalrRNA.gtf
- Sorted bam files in case that helped (tried both sorted and unsorted but didn't seem to matter)
samtools sort tophat_out_Indx26_RanHex/accepted_hits.bam tophat_out_Indx26_RanHex.sorted.bam samtools sort tophat_out_Indx27_dT/accepted_hits.bam tophat_out_Indx27_dT.sorted samtools sort tophat_out_Indx28_FISSEQRT/accepted_hits.bam tophat_out_Indx28_FISSEQRT.sorted samtools sort tophat_out_Indx29_Top48/accepted_hits.bam tophat_out_Indx29_Top48.sorted samtools sort tophat_out_Indx32_gDNA/accepted_hits.bam tophat_out_Indx32_gDNA.sorted
- Tried using intersectBed on bam file and gtf
- Tried sorted and unsorted bam file
- Was getting gibberish output until I figured out I needed -bed option to counter -abam option (which allows bam input but also makes bam output)
- Tried converting gtf to bed file (perl /home/mzcai/Perlscripts/gtf2bed.pl Homo_sapiens.GRCh37.75.totalrRNA.gtf > Homo_sapiens.GRCh37.75.totalrRNA.bed)
intersectBed -abam tophat_out_Indx26_RanHex.sorted.bam -b ~/Genomes/Homo_sapiens.GRCh37.75.totalrRNA.gtf > intersect_tophat_totalrRNA_Indx26_RanHex.bed -bed
Result Summary
- No matter what I tried the output file is always blank/empty, even if I use full gtf file (~/Genomes/Homo_sapiens.GRCh37.75.gtf)
- Must be something wrong with my command or intersectBed because there should definitely be overlap with the full genome annotation gtf