Editing
Matt:LabNotes/2014-7-14
Jump to navigation
Jump to search
Warning:
You are not logged in. Your IP address will be publicly visible if you make any edits. If you
log in
or
create an account
, your edits will be attributed to your username, along with other benefits.
Anti-spam check. Do
not
fill this in!
==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 **[[Media:Seq_Tophat_rRNA_Analysis.xlsx | 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 =====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: {| {{table}} | align="center" style="background:#f0f0f0;"|'''''' | align="center" style="background:#f0f0f0;"|'''Total aligned''' | align="center" style="background:#f0f0f0;"|'''Annotated with ignore.strand=FALSE''' | align="center" style="background:#f0f0f0;"|'''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 |} ====Bioconductor with TranscriptDb from Biomart==== *Same as Bioconductor with TranscriptDb from UCSC except get TranscriptDb from Biomart with this line: txdb=makeTranscriptDbFromBiomart(biomart='ensembl',dataset='hsapiens_gene_ensembl',transcript_ids=NULL,circ_seqs=DEFAULT_CIRC_SEQS,filters="",id_prefix="ensembl_",host="www.biomart.org",port=80,miRBaseBuild=NA) =====Results===== *Got 0 hits for every gene **Tried Indx32_gDNA and Indx26_RanHex ====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 *Continued on: [[Matt:LabNotes/2014-7-24]]
Summary:
Please note that all contributions to ZhangLabWiki may be edited, altered, or removed by other contributors. If you do not want your writing to be edited mercilessly, then do not submit it here.
You are also promising us that you wrote this yourself, or copied it from a public domain or similar free resource (see
ZhangLabWiki:Copyrights
for details).
Do not submit copyrighted work without permission!
Cancel
Editing help
(opens in new window)
Template used on this page:
Template:Table
(
edit
)
Navigation menu
Personal tools
Not logged in
Talk
Contributions
Create account
Log in
Namespaces
Page
Discussion
English
Views
Read
Edit
View history
More
Search
Navigation
Main Page
Current events
Recent changes
Random page
Investigators
Matt Cai
Song Chen
Eric Chu
Dinh Diep
Elizabeth Duong
Shicheng Guo
Alan Fung
Daniel Jacobsen
Blue Lake
Huy Lam
Alice Li
Andrew Richards
Brandon Sos
Chris Wei
Yan Wu
Kun Zhang
Tools
What links here
Related changes
Special pages
Page information