Rui:Data Analysis: Difference between revisions
>RuiLiu (→Tophat) |
>RuiLiu |
||
Line 576: | Line 576: | ||
* samtools view -h accepted_hits.bam > accepted_hits.sam | * 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 & | * 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 | |||
* Why no NGN3???? |
Revision as of 18:58, 9 August 2011
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]
Haplotyping 6.21.11
Data
ruiliu@genome-miner:~$ cd /home/kunzhang/haplotyping/Data/HL095/
ruiliu@genome-miner:/home/kunzhang/haplotyping/Data/HL095$ ls -l s_4_Indx1.*
-rw-r--r-- 1 kunzhang kunzhang 114237942 2011-06-13 22:07 s_4_Indx1.bowtie.pileup -rw-r--r-- 1 kunzhang kunzhang 39977658 2011-06-13 21:28 s_4_Indx1.bowtie.sorted.bam -rw-r--r-- 1 kunzhang kunzhang 4442210 2011-06-13 22:05 s_4_Indx1.bowtie.sorted.unique.bam -rw-r--r-- 1 kunzhang kunzhang 278 2011-06-13 21:29 s_4_Indx1.bowtie.variants.txt -rw-r--r-- 1 kunzhang kunzhang 7103 2011-06-14 15:58 s_4_Indx1.hapCompare.txt -rw-r--r-- 1 kunzhang kunzhang 9399 2011-06-14 16:56 s_4_Indx1.idioGraph.txt -rw-r--r-- 1 kunzhang kunzhang 675115 2011-06-13 21:29 s_4_Indx1.log -rw-r--r-- 1 kunzhang kunzhang 273522562 2011-06-13 18:08 s_4_Indx1.txt
ruiliu@genome-miner:/home/kunzhang/haplotyping/Data/HL095$ less s_4_Indx1.bowtie.pileup
1 557399 C A 25 25 25 1 ^:A f 1 557400 A A 25 0 25 1 . f 1 557401 T T 25 0 25 1 . f 1 557402 C C 25 0 25 1 . f 1 557403 A A 25 0 25 1 . c 1 557404 C C 25 0 25 1 . g 1 557405 T T 25 0 25 1 . g 1 557406 A A 25 0 25 1 . g 1 557407 G G 25 0 25 1 . g 1 557408 A A 25 0 25 1 . e 1 557409 C C 25 0 25 1 . f
ruiliu@genome-miner:/home/kunzhang/haplotyping/Data/HL095$ more s_4_Indx11.log
/home/kunzhang/softwares/bowtie-latest/bowtie -k 1 -l 32 -n 2 -m 1 --best --quiet -p 4 /home/kunzhang/HsGenome/1KG.ref/HsB36m s_4_Indx 11.txt s_4_Indx11.bowtie.map.1 /home/kunzhang/Exome/Data/scripts/trimTerminalMismatches.pl < s_4_Indx11.bowtie.map.1 > s_4_Indx11.bowtie.map.2 /home/kunzhang/softwares/samtools-0.1.7_x86_64-linux/bowtie2sam.pl s_4_Indx11.bowtie.map.2 > s_4_Indx11.bowtie.sam /home/kunzhang/softwares/samtools-0.1.7_x86_64-linux/samtools import /home/kunzhang/HsGenome/1KG.ref/human_b36_male.fa.gz.fai s_4_Indx 11.bowtie.sam s_4_Indx11.bowtie.bam [sam_header_read2] 114 sequences loaded. /home/kunzhang/softwares/samtools-0.1.7_x86_64-linux/samtools sort s_4_Indx11.bowtie.bam s_4_Indx11.bowtie.sorted java -Xmx2g -jar /home/kunzhang/softwares/picard-tools-1.38/MarkDuplicates.jar ASSUME_SORTED=TRUE REMOVE_DUPLICATES=TRUE INPUT=s_4_In dx11.bowtie.sorted.bam OUTPUT=s_4_Indx11.bowtie.sorted.unique.bam METRICS_FILE=dup.metric [Tue Jun 14 15:23:03 PDT 2011] net.sf.picard.sam.MarkDuplicates INPUT=s_4_Indx11.bowtie.sorted.bam OUTPUT=s_4_Indx11.bowtie.sorted.uni que.bam METRICS_FILE=dup.metric REMOVE_DUPLICATES=true ASSUME_SORTED=true MAX_SEQUENCES_FOR_DISK_READ_ENDS_MAP=50000 MAX_FILE_HANDL ES_FOR_READ_ENDS_MAP=8000 READ_NAME_REGEX=[a-zA-Z0-9]+:[0-9]:([0-9]+):([0-9]+):([0-9]+).* OPTICAL_DUPLICATE_PIXEL_DISTANCE=100 TMP_DIR =/tmp/kunzhang VERBOSITY=INFO QUIET=false VALIDATION_STRINGENCY=STRICT COMPRESSION_LEVEL=5 MAX_RECORDS_IN_RAM=500000 CREATE_INDEX=fals e CREATE_MD5_FILE=false INFO 2011-06-14 15:23:03 MarkDuplicates Start of doWork freeMemory: 374846832; totalMemory: 377225216; maxMemory: 1908932608 INFO 2011-06-14 15:23:03 MarkDuplicates Reading input file and constructing read end information. INFO 2011-06-14 15:23:03 MarkDuplicates Will retain up to 7575129 data points before spilling to disk. INFO 2011-06-14 15:23:03 MarkDuplicates Assuming input is coordinate sorted. INFO 2011-06-14 15:23:06 MarkDuplicates Read 887559 records. 0 pairs never matched. INFO 2011-06-14 15:23:07 MarkDuplicates After buildSortedReadEndLists freeMemory: 733397640; totalMemory: 925368320; maxMemory : 1908932608 INFO 2011-06-14 15:23:07 MarkDuplicates Will retain up to 59654144 duplicate indices before spilling to disk. INFO 2011-06-14 15:23:07 MarkDuplicates Traversing read pair information and detecting duplicates. INFO 2011-06-14 15:23:07 MarkDuplicates Traversing fragment information and detecting duplicates. INFO 2011-06-14 15:23:09 MarkDuplicates Sorting list of duplicate records. INFO 2011-06-14 15:23:09 MarkDuplicates After generateDuplicateIndexes freeMemory: 923987040; totalMemory: 1410072576; maxMemo ry: 1908932608 INFO 2011-06-14 15:23:09 MarkDuplicates Marking 837865 records as duplicates. INFO 2011-06-14 15:23:09 MarkDuplicates Found 0 optical duplicate clusters. INFO 2011-06-14 15:23:12 MarkDuplicates Before output close freeMemory: 920285040; totalMemory: 1408237568; maxMemory: 1908932 608 INFO 2011-06-14 15:23:12 MarkDuplicates After output close freeMemory: 876964080; totalMemory: 1364918272; maxMemory: 19089326 08 [Tue Jun 14 15:23:12 PDT 2011] net.sf.picard.sam.MarkDuplicates done. Runtime.totalMemory()=1364918272 /home/kunzhang/softwares/samtools-0.1.7_x86_64-linux/samtools pileup -c -f /home/kunzhang/HsGenome/1KG.ref/human_b36_male.fa s_4_Indx1 1.bowtie.sorted.unique.bam > s_4_Indx11.bowtie.pileup Use of uninitialized value in print at /home/kunzhang/Exome/Data/scripts/pileup2variants.pl line 145. Use of uninitialized value in print at /home/kunzhang/Exome/Data/scripts/pileup2variants.pl line 145.
ruiliu@genome-miner:/home/kunzhang/haplotyping/Data/HL095$ /home/kunzhang/softwares/samtools-0.1.7_x86_64-linux/samtools view s_4_Indx11.bowtie.sorted.bam | wc -l
887559
ruiliu@genome-miner:/home/kunzhang/haplotyping/Data/HL095$ /home/kunzhang/softwares/samtools-0.1.7_x86_64-linux/samtools view s_4_Indx11.bowtie.sorted.unique.bam | wc -l
49694
Reads
Barcode | Reads | Unique | Percentage |
Indx1 | 163556 | 15899 | 9.7 |
Indx2 | 196625 | 16718 | 8.5 |
Indx3 | 127140 | 10123 | 8.0 |
Indx4 | 189464 | 13511 | 7.1 |
Indx5 | 9536 | 4856 | 50.9 |
Indx6 | 220670 | 12779 | 5.8 |
Indx7 | 185404 | 12479 | 6.7 |
Indx8 | 209445 | 14185 | 6.8 |
Indx9 | 175999 | 12100 | 6.9 |
Indx10 | 212683 | 13919 | 6.5 |
Indx11 | 188707 | 14591 | 7.7 |
Indx12 | 326082 | 16300 | 5.0 |
RNAseq analysis from 7.18.11
- 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
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
Bowtie Index
- 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
- 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
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
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
Summary on Tophat/cufflinks/cuffcompare/cuffdiff
Tophat
- 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
- 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% |
RNAseq analysis on HL098
Database download
iGenome collection [4], save in ~/iGenome/
- wget ftp://igenome:G3nom3s4u@ftp.illumina.com/Homo_sapiens/UCSC/hg19/Homo_sapiens_UCSC_hg19.tar.gz
- wget ftp://igenome:G3nom3s4u@ftp.illumina.com/Mus_musculus/UCSC/mm9/Mus_musculus_UCSC_mm9.tar.gz
- tar -zxvf Homo_sapiens_UCSC_hg19.tar.gz
- tar -zxvf Mus_musculus_UCSC_mm9.tar.gz
- 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
- Why no NGN3????