Rui:RNAseq analysis from 7.18.11
Jump to navigation
Jump to search
RNAseq analysis from 7.18.11[edit]
- ID1 and ID2 are human samples; ID3-6 are mouse samples. De-multiplexing will give clearer idea on library mapping.
- Combine them into a single file, I can first check the sequencing quality with Galaxy.
- Athurva converted qseq files (55bp) for lane 1 (120) into a single fastq file [ruiliu@genome-miner:~/RNAseq/s_1_1_sequence.txt]
- s_1_1_sequence.txt: 4.4 GB, 50 million lines [wc -l s_1_1_sequence.txt]
- [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)
- 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??)
- 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:
- Which one is the mappable reads?
- 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?
- 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?.
- 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]
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% |