Matt:LabNotes/2014-7-14: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Mzcai
>Mzcai
Line 113: Line 113:
| gDNA||44982||15441||31096
| 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====
====Bedtools intersect====

Revision as of 22:10, 16 July 2014

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:
    • 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
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

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