Rui: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]
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 | ' | ' | ' | bowtie | ' | ' | ' | ' | removeClonalHits.pl | ' | ' | ' | ' | ' | ' | ' |
Sample | file size | fastq reads in | reads out | Genome | Accepted hits | 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% |