Rui:Data Analysis

From ZhangLabWiki
Revision as of 23:15, 14 October 2011 by >RuiLiu (→‎Data analysis)
Jump to navigation Jump to search

Examples

Jeff's samples:
Kun:LabNotes/Haplotyping/2010-9-17 [1]
Transcriptome:
Kun:LabNotes/SingleCellExpr/2011-6-13 [2]
Haplotyping:
Kun:LabNotes/Haplotyping/2011-6-14 [3]

Data analysis

RNAseq analysis on HL098

Database download

iGenome collection [4], save in ~/iGenome/

  • BowtieIndex: ~/iGenome/Homo_sapiens/UCSC/hg19/Sequence/BowtieIndex | ~/iGenome/Mus_musculus/UCSC/mm9/Sequence/BowtieIndex
  • Annotation: ~/iGenome/Homo_sapiens/UCSC/hg19/Annotation/Genes/genes.gtf

Tophat

nohup tophat -p 6 --solexa1.3-quals -o ./tophat-g/ -r 250 --mate-std-dev 80 -G ~/iGenome/Homo_sapiens/UCSC/hg19/Annotation/Genes/genes.gtf ~/iGenome/Homo_sapiens/UCSC/hg19/Sequence/BowtieIndex/genome ./seq/s_1_1_Indx1.txt ./seq/s_1_2_Indx1.mod.txt

nohup tophat -p 6 --solexa1.3-quals -o ./tophat-g/ -r 250 --mate-std-dev 80 -G ~/iGenome/Homo_sapiens/UCSC/hg19/Annotation/Genes/genes.gtf ~/iGenome/Homo_sapiens/UCSC/hg19/Sequence/BowtieIndex/genome ./seq/s_1_1_Indx2.txt ./seq/s_1_2_Indx2.mod.txt

Samtools and clonalreads removal

samtools flagstat accepted_hits.bam

~/RNAtools/removalclonalreads.pl accepted_hits.bam


Pair-end mapping ' ' ' ' ' Tophat ' ' ' ' ' removeClonalHits.pl ' ' ' ' ' ' '
file size fastq reads in reads out % options accepted hits reads hits pair/single Properly paired % of total mapped reads % of total pro-rm clonal reads % of total % of mappable reads then uniquely mapped % of total potential coverage
s_1_Indx1 s_1_1_Indx1 left read 2,307,029 2,300,738 0.9973 r:250 1,705,993 734,602 1,122,540 0.657998011 1,615,384 0.3501 1,191,029 0.2581 0.7373 1,160,079 0.2514 0.04640316
s_1_2_Indx1 right read 2,307,029 2,304,795 0.9990 STD:80 971,391 583,453 941,458 0.551853378
s_1_Indx1 s_1_1_Indx1 left read 2,307,029 2,300,738 0.9973 (+g) 1,733,904 746,989 1,156,628 0.667065766 1,643,007 0.3561 1,218,033 0.2640 0.7413 1,186,899 0.2572 0.04747596
s_1_2_Indx1 right read 2,307,029 2,304,795 0.9990 UCSC 986,915 577,276 965,638 0.556915492
s_1_Indx2 s_1_1_Indx2 left read 3,078,335 3,069,878 0.9973 r:250 3,882,019 1,679,990 2,750,090 0.70841745 3,678,987 0.5976 2,642,026 0.4291 0.7181 2,585,255 0.4199 0.1034102
s_1_2_Indx2 right read 3,078,335 3,075,318 0.9990 STD:80 2,202,029 1,131,929 2,293,768 0.590869854
s_1_Indx2 s_1_1_Indx2 left read 3,078,335 3,069,878 0.9973 (+g) 3,913,815 1,694,018 2,792,224 0.713427691 3,711,732 0.6029 2,674,065 0.4343 0.7204 2,617,136 0.4251 0.10468544
s_1_2_Indx2 right read 3,078,335 3,075,318 0.9990 UCSC 2,219,797 1,121,591 2,344,834 0.59911723

Cufflinks and cuffdiff

samtools view -h ../accepted_hits.bam > accepted_hits.sam (have to convert to sam first, i can't feed cufflinks with bam file --- error: segment fault...)

cufflinks -g ~/iGenome/Homo_sapiens/UCSC/hg19/Annotation/Genes/genes.gtf -m 250 -s 80 accepted_hits.sam (-g option is probably unnecessary...)


RNAseq analysis on HL099

Download files from SeqStore

  • transfer lane 2 to 6 files from HL099_2 directory in SeqStore/genome-minor
  • cat s_*_Indx1.txt > Indx1.txt
  • cat s_*_Indx2.txt > Indx2.txt
Libraries (Indx1-6) in HL098 lane 1 and HL099 lane 2 are same. Libraries in HL099 lane 3 to 6 are same, re-adjusted ratio based on Indx2 (set as 1)
ID Sample file size relative ratio
HL098 s_1_Indx1.txt 1,063,078,988 0.75
HL098 s_1_Indx2.txt 1,418,503,096 1.00
HL098 s_1_Indx3.txt 1,126,757,636 0.79
HL098 s_1_Indx4.txt 1,773,739,660 1.25
HL098 s_1_Indx5.txt 4,002,812,404 2.82
HL098 s_1_Indx6.txt 4,927,760,680 3.47
10.09
HL099 s_2_Indx1.txt 799,803,430 0.73
HL099 s_2_Indx2.txt 1,088,795,544 1.00
HL099 s_2_Indx3.txt 816,515,216 0.75
HL099 s_2_Indx4.txt 668,702,922 0.61
HL099 s_2_Indx5.txt 1,263,996,032 1.16
HL099 s_2_Indx6.txt 1,152,064,234 1.06
5.32
HL099 s_3_Indx1.txt 2,055,529,230 1.58
HL099 s_3_Indx2.txt 1,300,750,376 1.00
HL099 s_3_Indx3.txt 1,054,455,616 0.81
HL099 s_3_Indx4.txt 1,205,324,668 0.93
HL099 s_4_Indx1.txt 2,075,122,264 1.59
HL099 s_4_Indx2.txt 1,308,387,872 1.00
HL099 s_4_Indx3.txt 1,057,664,912 0.81
HL099 s_4_Indx4.txt 1,212,653,382 0.93
HL099 s_5_Indx1.txt 1,942,428,848 1.58
HL099 s_5_Indx2.txt 1,225,815,974 1.00
HL099 s_5_Indx3.txt 984,527,600 0.80
HL099 s_5_Indx4.txt 1,129,401,814 0.92
HL099 s_6_Indx1.txt 2,046,395,112 1.58
HL099 s_6_Indx2.txt 1,292,688,682 1.00
HL099 s_6_Indx3.txt 1,044,933,252 0.81
HL099 s_6_Indx4.txt 1,189,162,084 0.92


Tophat

  • nohup tophat -p 6 --solexa1.3-quals -G ~/iGenome/Homo_sapiens/UCSC/hg19/Annotation/Genes/genes.gtf ~/iGenome/Homo_sapiens/UCSC/hg19/Sequence/BowtieIndex/genome ./seq/Indx1.txt
  • nohup tophat -p 6 --solexa1.3-quals -G ~/iGenome/Homo_sapiens/UCSC/hg19/Annotation/Genes/genes.gtf ~/iGenome/Homo_sapiens/UCSC/hg19/Sequence/BowtieIndex/genome ./seq/Indx2.txt
  • I also tried to feed tophat with the junctions file generated in HL098 (w/ PE mapping strategy) to improve accuracy. However, the bed_to_juncs < junctions.bed > new_list.juncs command didn't work at all! (http://tophat.cbcb.umd.edu/manual.html)
  • samtools flagstat accepted_hits.bam --> count the total hits
  • ~/RNAtools/removalClonalhits.pl accepted_hits.bam --> total reads based on hits, reads after clonal reads removal, and reads after clonal reads removal plus uniquely mapped in genome (considered as mappable reads)
  • potential coverage is based on estimation in the ENCODE standard: mappable reads (uniquely mapped un-clonal reads)/ (2X10^9 mRNAs/cell X ~1500bp/mRNA)
  • Big problem in this run (HL099) is clonal reads -- more than 50% of total mapped reads! (the one in HL098 is roughly 28-30%)


run Indx file size QC ' ' Tophat removeClonalHits.pl ' ' ' ' ' ' '
reads-in reads-out % accepted hits mapped reads % of reads-in pro-rm clonal reads % of reads-in % of mapped reads then uniquely mapped % of reads-in potential coverage
HL099 Indx1 8,919,278,884 46,366,051 46,356,341 0.9998 23,635,554 20,388,632 0.44 8,637,058 0.19 0.42 8,355,371 0.18 0.33421484
HL099 Indx2 6,216,438,448 32,314,674 32,306,989 0.9998 27,114,017 23,424,639 0.72 10,904,543 0.34 0.47 10,567,772 0.33 0.42271088

Cuffdiff

  • I never get the original command (listed in cufflinks website) worked! Just as when I run cufflinks, I have to convert bam to sam before running, otherwise, it will show segmentation fault.
The original command
  • Run Cuffdiff: Take the annotated transcripts for your genome (as GFF or GTF) and provide them to cuffdiff along with the BAM files from TopHat for each replicate: cuffdiff annotation.gtf mock_rep1.bam,mock_rep2.bam \ knockdown_rep1.bam,knockdown_rep2.bam
Run with accepted_hits.bam from tophat
  • nohup cuffdiff ~/iGenome/Homo_sapiens/UCSC/hg19/Annotation/Genes/genes.gtf ../Indx1/tophat_out/accepted_hits.bam ../Indx2/tophat_out/accepted_hits.bam &
  • [1]+ Segmentation fault nohup cuffdiff ~/iGenome/Homo_sapiens/UCSC/hg19/Annotation/Genes/genes.gtf ../Indx1/tophat_out/accepted_hits.bam ../Indx2/tophat_out/accepted_hits.bam
Run with sorted.acceptedHits.sam from removalClonalhits.pl
  • nohup cuffdiff ~/iGenome/Homo_sapiens/UCSC/hg19/Annotation/Genes/genes.gtf ../Indx1/tophat_out/sorted.acceptedHits.sam ../Indx2/tophat_out/sorted.acceptedHits.sam &
  • Error: this SAM file doesn't appear to be correctly sorted!
       current hit is at chrX:64289, last one was at chrM:16509
Cufflinks requires that if your file has SQ records in
the SAM header that they appear in the same order as the chromosomes names 
in the alignments.
If there are no SQ records in the header, or if the header is missing,
the alignments must be sorted lexicographically by chromsome
name and by position.
Convert bam to sam with header, then cuffdiff
  • samtools view -h accepted_hits.bam > accepted_hits.sam
  • nohup cuffdiff ~/iGenome/Homo_sapiens/UCSC/hg19/Annotation/Genes/genes.gtf ~/HL099/Indx1/tophat_out/cuffcompare/accepted_hits.sam ~/HL099/Indx2/tophat_out/cuffcompare/accepted_hits.sam &
  • gene_exp.diff file has ~1700 genes with significant... repeat direct-cuffdiff on HL098 to compare [636 genes with significance as same as ones using tophat-cufflinks-cuffcompare-cuffdiff methods]
  • Why no NGN3????