Rui:RNAseq analysis from 7.18.11

From ZhangLabWiki
Jump to navigation Jump to search

RNAseq analysis from 7.18.11[edit]

  1. ID1 and ID2 are human samples; ID3-6 are mouse samples. De-multiplexing will give clearer idea on library mapping.
  2. Combine them into a single file, I can first check the sequencing quality with Galaxy.
  3. Athurva converted qseq files (55bp) for lane 1 (120) into a single fastq file [ruiliu@genome-miner:~/RNAseq/s_1_1_sequence.txt]
  4. s_1_1_sequence.txt: 4.4 GB, 50 million lines [wc -l s_1_1_sequence.txt]
  5. [head -n 1000000 s_1_1_sequence.txt] for the firstmillion lines [ruiliu@genome-miner:~/RNAseq/firstmillion.txt] uploaded to Galaxy (up limit is 2GB)
  6. Groomer --> fastq summary statistics (250000 fastq reads) --> compute quality statistics (good) --> draw quality score boxplot (40-35) --> Tophat (accepted hits /splice junction) --> flagstat (143296 in total, 0 QC failure, 0 duplicates, 143296 mapped (100.00%), 0 paired in sequencing, 0 read1) --> Mark Duplicate reads (empty??)
  7. To upload the original file (4GB), zip it to [tar -zcvf ~/new.tar.gz s_1_sequence.txt] new.tar.gz under home directly (permission issue), download to my laptop and install Cyberduck on Mac to upload file via FTP [main.g2.bx.psu.edu] with username [bsos@ucsd.edu] and password [zhanglab]

overview of reads[edit]

Questions:

  1. Which one is the mappable reads?
  2. Clonal reads can't be removed if the expression level is compared, right? So stopping PCR at exponential stage is critical. Low expressed RNAs will be sacrificed to high expressed RNAs if over-amplifying RNAs and total input amount is limited for sequencing. And eventually, increasing clonal reads will decrease sequencing coverage, right? But still, PCR amplification plus clustering could also lead inaccuracy in expression level estimate, right?
  3. Percentage of reads after clonal reads removal to mappable reads in Ind1 to Ind4 is similar (76%; 76.7%; 76.8%;78%), which means these libraries are almost evenly amplified, right?.
  4. Multiple loci reads also can't be removed b/c short reads of coding sequence are possibly mapped to paralogs, right?
110714_HL098 ' ' FASTQ stat. ' ' TopHat ' ' removeClonalHits.pl ' ' ' ' '
ID Sample file size fastq reads ASCII range Decimal range Genome Accepted hits Junction mapped reads % of total pro-rm clonal reads % of total then uniquely mapped % of total
s_1_Indx1 GFP+ 531.5MB 2,307,029 #\'(35) - \'I\'(73) 2-40 Hs 19 889,766 36,435 718,472 0.3114 546,308 0.236801531 511,691 0.221796518
s_1_Indx2 GFP- 709.3MB 3,078,335 \'B\'(66) - \'h\'(104) 33-71 Hs 19 1,963,553 61,496 1,615,433 0.5248 1,237,498 0.402002381 1,170,098 0.380107428
s_1_Indx3 E9.5 563.4MB 2,445,243 \'B\'(66) - \'h\'(104) 33-71 mm 9 1,983,651 47,929 1,491,621 0.6100 1,143,810 0.467769461 1,039,861 0.425258758
s_1_Indx4 E11.5 886.9MB 3,849,227 \'B\'(66) - \'h\'(104) 33-71 mm 9 1,618,640 44,625 1,208,317 0.3139 951,936 0.247305758 856,563 0.222528575


Run Tophat in meangenemachine[edit]

Bowtie Index[edit]

  • Download Hs19
 Wget ftp://ftp.ncbi.nih.gov/genomes/H_Sapine/Assambled_chromosome/hs_ref_GRCh37.p2*.fa.gz
  • Combine all together
 less *.gz > hs_ref_GRCh37.fa
  • Bowtie Index
 bowtie-build hs_ref_GRCh37.fa hs_ref_GRCh37 > log.blablabla

Sequence files[edit]

  • check fastq file with ASCII table to see if 33 code or 64 code
  • combine 2 ends files into 1 file (somehow, less command is not working well in Mac)
 cat file1.gz file2.gz | gunzip > new file
  • substitute "#0/3" with "#03" (reason: position confusion in the sorted.acceptedhit.sam file generated from TopHat,esp. running removalclonalread.ps after that; more information for sed : http://www.grymoire.com/Unix/Sed.html#uh-0)
 sed 's:#0/3:#03:g' file > newfile
 sed 's/#0\/3/#03/g' file > newfile
  • run QC for each file and download
 ~/RNAseqTools/fastqc file
  • run Tophat
 in /BowtieIndices/ to find the human or mm ref
 tophat -p 4 --solexa1.3-quals /media/1TB_store1/BowtieIndices/mm_ref_MGSCv37 file &
 For pair-end:
 tophat -p 4 --solexa1.3-quals -r "distance" --mate-std-dev "STDEV" -o "output" /media/1TB_store1/BowtieIndices/hs_ref_GRCh37.p2 file &

Tophat summary[edit]

110714_HL098 ' ' TopHat ' ' ' ' ' removeClonalHits.pl ' ' ' ' ' '
ID Sample file size fastq reads in reads out Genome Accepted hits Junction mapped reads % of total pro-rm clonal reads % of total % of mappable reads then uniquely mapped % of total
s_1_Indx1 GFP+ 1,063,078,988 4,614,058 4,605,533 0.9982 Hs 19 1,819,753 1,646,034 0.3567 1,191,550 0.2582 0.723891487 1,140,325 0.2471
s_1_Indx2 GFP- 1,418,503,096 6,156,670 6,145,196 0.9981 Hs 19 4,156,441 3,742,990 0.6080 2,637,547 0.428404803 0.704663117 2,543,973 0.413206003
s_1_Indx3 E9.5 1,126,757,636 4,890,486 4,881,588 0.9982 mm 9 4,483,009 3,529,559 0.7217 2,350,598 0.480647118 0.665974984 2,192,476 0.448314544
s_1_Indx4 E11.5 1,773,739,660 7,698,454 7,685,366 0.9983 mm 9 3,537,587 2,757,814 0.3582 1,915,170 0.248773325 0.6944522 1,771,219 0.230074636
s_1_Indx5 E13.5m 4,002,812,404 17,373,520 17,342,655 0.9982 mm 9 6,601,730 5,227,070 0.3009 3,214,319 0.1850 0.614937049 3,023,072 0.1740
s_1_Indx6 E13.5f 4,927,760,680 21,387,852 21,349,698 0.9982 mm 9 5,999,758 4,744,699 0.221840837 2,967,058 0.13872632 0.625341671 2,783,215 0.130130646


Cell numbers[edit]

ID Sample RNA amount Cell numbers % of mappable read to total
s_1_Indx1 GFP+ 20ng 0.2471
s_1_Indx2 GFP- 12ng 0.413206003
s_1_Indx3 E9.5 ~316 cells 0.448314544
s_1_Indx4 E11.5 ~2000 cells 0.230074636
s_1_Indx5 E13.5m ~3000 cells 0.1740
s_1_Indx6 E13.5f ~3000 cells 0.130130646

Flowchart[edit]

File:Flowchart 1.jpg

File:Flowchart 2.jpg

Summary on Tophat/cufflinks/cuffcompare/cuffdiff[edit]

Tophat[edit]

  • reads results from single end mapping and from paired end mapping are quite similar; parameters such as internal-length and STDEV have least, if not no, effect.
  • Paired end mapping provide more accurate information, and show significant difference when using cufflinks
' ' ' TopHat ' ' ' ' ' ' ' ' removeClonalHits.pl ' ' ' ' ' ' '
Sample file size fastq reads in reads out Genome Accepted hits hits from each read mapped pairs (up) mapped single (down) properly located pairs % 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 GFP+ 1,063,078,988 4,614,058 4,605,533 0.9982 Hs 19 1,819,753 1,646,034 0.3567 1,191,550 0.2582 0.7239 1,140,325 0.2471 0.045613
s_1_Indx2 GFP- 1,418,503,096 6,156,670 6,145,196 0.9981 Hs 19 4,156,441 3,742,990 0.6080 2,637,547 0.4284 0.7047 2,543,973 0.4132 0.10175892
s_1_Indx1 s_1_1_Indx1 left read 2,307,029 2,300,738 0.9973 r:225 1,706,112 734,674 1,122,672 0.658029485 1,615,377 0.3501 1,191,056 0.2581 0.7373 1,160,088 0.2514 0.04640352
s_1_2_Indx1 right read 2,307,029 2,304,795 0.9990 STD:212 971,438 583,440 955,602 0.560105081
s_1_Indx2 s_1_1_Indx2 left read 3,078,335 3,069,878 0.9973 3,881,050 1,679,506 2,749,076 0.708333054 3,678,975 0.5976 2,642,078 0.4291 0.7182 2,585,404 0.4199 0.10341616
s_1_2_Indx2 right read 3,078,335 3,075,318 0.9990 2,201,544 1,131,974 2,346,216 0.604531248
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_Indx2 s_1_1_Indx2 left read 3,078,335 3,069,878 0.9973 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 2,202,029 1,131,929 2,293,768 0.590869854

Cufflinks[edit]

  • cufflinks: processing loci greatly differ in SE/PE mapping and correction with GTF (tip: RefFlast from Refseq in UCSC, can't used gene_id from Ensemble, if gene name is expected in output files)
  • Correction wit GTF largely reduced loci input and # of output in genes/isoforms (cufflinks, cuffcompare and cuffdiff), probably due to better annotation? (guess, can't find details in manual or paper...)
  • No much parameter options for cuffcompare or cuffdiff
  • Odd: different gene_exp in tophat_225_212 and tophat_250_80 are exactly same, although they are different in any other results.
' ' cufflinks ' ' ' cuffcompare ' ' ' ' ' cuffdiff '
loci total map mass wc -l Missed exons Wrong exons Missed introns Wrong introns Missed loci Wrong loci gene isoform
tophat_SE Indx1 131,431 1645874.66 66,396 42.60% 7.40% 48.50% 0.30% 0.00% 23.40% 1745 4623
Indx2 203,286 3742794.31 117,273 37.90% 8.50% 42.80% 0.60% 0.00% 27.10%
tophat_PE Indx1_left 95,470 1076008.63 29,623 43.20% 5.60% 48.40% 0.20% 0.00% 19.90% 332* 646
225/212 Indx1_right (-g) 40,331 441,677
Indx2_left 39.50% 5.30% 43.80% 0.20% 0.00% 19.90%
Indx2_right
tophat_PE Indx1_left
250/80 Indx1_right (-g) 40,341 1076014.17 441,831 0.00% 0.30% 0.20% 0.00% 0.00% 2.50% 332* 1462
Indx2_left
Indx2_right (-g) 63,059 2361671.63 459,770 0.00% 0.60% 0.20% 0.10% 0.00% 4.10%