Rui:Data Analysis: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>RuiLiu
>RuiLiu
 
(35 intermediate revisions by the same user not shown)
Line 1: Line 1:
==Examples==
==Examples==
KUN's pipeline:
http://genome-tech.ucsd.edu/LabNotes/index.php/RNAseq_pipeline [http://genome-tech.ucsd.edu/LabNotes/index.php/RNAseq_pipeline]
  Jeff's samples:
  Jeff's samples:
  Kun:LabNotes/Haplotyping/2010-9-17 [http://genome-tech.ucsd.edu/LabNotes/index.php/Kun:LabNotes/Haplotyping/2010-9-17]
  Kun:LabNotes/Haplotyping/2010-9-17 [http://genome-tech.ucsd.edu/LabNotes/index.php/Kun:LabNotes/Haplotyping/2010-9-17]
Line 9: Line 12:
  Kun:LabNotes/Haplotyping/2011-6-14 [http://genome-tech.ucsd.edu/LabNotes/index.php/Kun:LabNotes/Haplotyping/2011-6-14]
  Kun:LabNotes/Haplotyping/2011-6-14 [http://genome-tech.ucsd.edu/LabNotes/index.php/Kun:LabNotes/Haplotyping/2011-6-14]


==Haplotyping 6.21.11==
==Data analysis==
 
*[[Rui:Haplotyping 6.21.11|Haplotyping 6.21.11]]
===Data===
*[[Rui:RNAseq analysis from 7.18.11|RNAseq analysis from 7.18.11]]
ruiliu@genome-miner:~$ cd /home/kunzhang/haplotyping/Data/HL095/
*[[Rui:RNAseq analysis on HL098|RNAseq analysis on HL098]]
 
*[[Rui:RNAseq analysis on HL099|RNAseq analysis on HL099]]
ruiliu@genome-miner:/home/kunzhang/haplotyping/Data/HL095$ ls -l s_4_Indx1.*
*[[Rui:RNAseq analysis on HL101|RNAseq analysis on HL101]]
-rw-r--r-- 1 kunzhang kunzhang 114237942 2011-06-13 22:07 s_4_Indx1.bowtie.pileup
*[[Rui:RNAseq analysis on Hiseq111005|RNAseq_PGCsamples on Hiseq111005]]
-rw-r--r-- 1 kunzhang kunzhang  39977658 2011-06-13 21:28 s_4_Indx1.bowtie.sorted.bam
*[[Rui:Hap analysis on HL108|Hap analysis on HL108]]
-rw-r--r-- 1 kunzhang kunzhang  4442210 2011-06-13 22:05 s_4_Indx1.bowtie.sorted.unique.bam
*[[Rui:RNAseq analysis on HL109|RNAseq analysis on HL109]]
-rw-r--r-- 1 kunzhang kunzhang      278 2011-06-13 21:29 s_4_Indx1.bowtie.variants.txt
*[[Rui:DNAseq analysis on HL114|DNAseq analysis on HL114]]
-rw-r--r-- 1 kunzhang kunzhang      7103 2011-06-14 15:58 s_4_Indx1.hapCompare.txt
*[[Rui:DNAseq analysis on Hiseq120313|DNAseq analysis on Hiseq120313]]
-rw-r--r-- 1 kunzhang kunzhang      9399 2011-06-14 16:56 s_4_Indx1.idioGraph.txt
*[[Rui:DNAseq analysis on HL118|DNAseq analysis on HL118]]
-rw-r--r-- 1 kunzhang kunzhang    675115 2011-06-13 21:29 s_4_Indx1.log
*[[Rui:totoRNAseq on 121226_HL140|totoRNAseq on 121226_HL140]]
-rw-r--r-- 1 kunzhang kunzhang 273522562 2011-06-13 18:08 s_4_Indx1.txt
*[[Rui:totoRNAseq on 130120_HL141|totoRNAseq on 130120_HL141]]
 
*[[Rui:totoRNAseq on 130220_HL144|totoRNAseq on 130220_HL144]]
ruiliu@genome-miner:/home/kunzhang/haplotyping/Data/HL095$ less s_4_Indx1.bowtie.pileup
*[[Rui:totoRNAseq on BL-UHRR-IVT-May28_2013|totoRNAseq on BL-UHRR-IVT-May28_2013]]
1      557399  C      A      25      25      25      1      ^:A    f
*[[Rui:totoRNAseq on 131009_HL161|totoRNAseq on 131009_HL161]]
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===
{| {{table}}
| align="center" style="background:#f0f0f0;"|'''Barcode'''
| align="center" style="background:#f0f0f0;"|'''Reads'''
| align="center" style="background:#f0f0f0;"|'''Unique'''
| align="center" style="background:#f0f0f0;"|'''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?
 
{| {{table}} border=1
| align="center" style="background:#f0f0f0;"|'''110714_HL098'''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|'''FASTQ stat.'''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|'''TopHat'''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|'''removeClonalHits.pl'''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
|-
| 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====
 
{| {{table}}
| align="center" style="background:#f0f0f0;"|'''110714_HL098'''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|'''TopHat'''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|'''removeClonalHits.pl'''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
|-
| 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====
{| {{table}} border=1
| align="center" style="background:#f0f0f0;"|'''ID'''
| align="center" style="background:#f0f0f0;"|'''Sample'''
| align="center" style="background:#f0f0f0;"|'''RNA amount'''
| align="center" style="background:#f0f0f0;"|'''Cell numbers'''
| align="center" style="background:#f0f0f0;"|'''% 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===
[[File:flowchart_1.jpg]]
 
[[File:flowchart_2.jpg]]
 
===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
 
{| {{table}} border=1
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|'''TopHat'''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|'''removeClonalHits.pl'''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
|-
| ||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.
 
{| {{table}} border=1
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|'''cufflinks'''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|'''cuffcompare'''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|'''cuffdiff'''
| align="center" style="background:#f0f0f0;"|''''''
|-
| ||||||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 [http://cufflinks.cbcb.umd.edu/igenomes.html], 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
 
 
{| {{table}} border=1
| align="center" style="background:#f0f0f0;"|'''Pair-end mapping'''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|'''Tophat'''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|'''removeClonalHits.pl'''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
|-
| ||||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)
 
{| {{table}} border=1
| align="center" style="background:#f0f0f0;"|'''ID'''
| align="center" style="background:#f0f0f0;"|'''Sample'''
| align="center" style="background:#f0f0f0;"|'''file size'''
| align="center" style="background:#f0f0f0;"|'''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===
==RNAseq analysis key notes==
* 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
* iGenome packages [http://tophat.cbcb.umd.edu/igenomes.html]
* 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
# UCSC genome package: [http://genome-tech.ucsd.edu/LabNotes/index.php/Rui:Tophat_mapping_10.27.11]
* samtools flagstat accepted_hits.bam --> count the total hits
# Ensemble genome package: [http://genome-tech.ucsd.edu/LabNotes/index.php/Rui:Tophat_11.21.11#Re-download_Ensemble_package_from_iGenome_.5B1.5D]
* ~/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)


{| {{table}} border=1
* Comparison of RNAseq analysis programs: [http://genome-tech.ucsd.edu/LabNotes/index.php/Rui:RNAseq_analysis_on_Hiseq111005#Comparison_of_program]
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|'''Tophat'''
| align="center" style="background:#f0f0f0;"|'''removeClonalHits.pl'''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
|-
| ||||||reads in||reads out||% ||accepted hits||mapped reads||% of total||pro-rm clonal reads||% of total||% of mappable reads||then uniquely mapped||% of total||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
|}

Latest revision as of 21:02, 31 October 2013

Examples[edit]

KUN's pipeline:
http://genome-tech.ucsd.edu/LabNotes/index.php/RNAseq_pipeline [1]
Jeff's samples:
Kun:LabNotes/Haplotyping/2010-9-17 [2]
Transcriptome:
Kun:LabNotes/SingleCellExpr/2011-6-13 [3]
Haplotyping:
Kun:LabNotes/Haplotyping/2011-6-14 [4]

Data analysis[edit]

RNAseq analysis key notes[edit]

  • iGenome packages [5]
  1. UCSC genome package: [6]
  2. Ensemble genome package: [7]
  • Comparison of RNAseq analysis programs: [8]