Rui:RNAseq analysis on Hiseq111005: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>RuiLiu
>RuiLiu
 
(32 intermediate revisions by the same user not shown)
Line 1: Line 1:
==RNAseq analysis on Hiseq111005==
==RNAseq analysis on Hiseq111005==
===Comparison of program===
{| {{table}} border=1
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|'''Meangenemachine'''
| align="center" style="background:#f0f0f0;"|'''Genome-miner'''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
|-
| Tophat||||v1.3.1||v1.3.1||||
|-
| Bowtie||||0.12.7.0||0.12.7.0||||
|-
| Samtools||||0.1.17||0.1.12a||||
|-
| BowtieIndex||||(UCSC)||(UCSC)||||(Ensemble)
|-
| ||.fa||2011.8.8||2011.8.6||2011.10.27||2011.11.21
|-
| ||.ebwt||2011.5.16||2011.5.16||2011.9.27||2011.9.27
|-
| Annotation||gene.gtf||2011.5.17||2011.5.17||2011.9.27||2011.9.27
|-
| ||file||~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/gene.gtf||/GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf||~/RNAseq/iGenome/UCSC/Mus_musculus/UCSC/mm9/Annotation/Genes/gene.gtf||[http://genome-tech.ucsd.edu/LabNotes/index.php/Rui:Tophat_11.21.11#Tophat_11.21.11]
|-
| cufflinks||||v1.0.3||v1.1.0 (2699)||||
|-
| cuffdiff||||v1.0.3||v1.1.0 (2699)||||
|}


===Data===
===Data===
Line 18: Line 48:
  scp ruiliu@genome-miner:~/RNAseq/Hiseq111005/Indx_seq/* ./
  scp ruiliu@genome-miner:~/RNAseq/Hiseq111005/Indx_seq/* ./


===Tophat w/o G correction===
===Tophat_10.25.11 w/o G correction===
*[[Rui:Tophat mapping|Tophat mapping]]
*[[Rui:Tophat mapping|Tophat mapping]]
*[[Rui:Transfer and convert files|Transfer and convert files]]
*[[Rui:Transfer and convert files|Transfer and convert files]]
*[[Rui:Cuffdiff|Cuffdiff]]
*[[Rui:Cuffdiff|Cuffdiff]]
*[[Rui:Cufflinks|Cufflinks]]
*[[Rui:Cufflinks_10.26.11|Cufflinks_10.26.11]]
*[[Rui:Matrix|Matrix]]
*[[Rui:Matrix|Matrix]]


===Tophat w/ G correction===
===Tophat_10.27.11 w/ G correction===
*[[Rui:Tophat mapping 10.27.11|Tophat mapping 10.27.11]]
*[[Rui:Tophat mapping 10.27.11|Tophat mapping 10.27.11]]
*[[Rui:Cufflinks 10.28.11|Cufflinks 10.28.11]]
*[[Rui:Cufflinks 10.28.11|Cufflinks 10.28.11]]
*[[Rui:Matrix 10.30.11|Matrix 10.30.11]]
*[[Rui:Matrix 10.30.11|Matrix 10.30.11]]
*[[Rui:Cufflinks_11.7.11|Cufflinks_11.7.11 w/ G annotation]]


==Repeat RNAseq analysis on Hiseq111005 from 11.10.11==
* Lesson 1: Use exactly same VERSION of programs/softwares to do analysis on ALL data in ONE computer
* Lesson 2: Random check for comparison
*[[Rui:Tophat_11.10.11|Tophat_11.10.11]]
*[[Rui:Cufflinks_11.11.11|Cufflinks_11.11.11]]
*[[Rui:Cuffdiff_11.11.11|Cuffdiff_11.11.11]]


====Tophat mapping====
* Tophat against UCSC.gtf
nohup tophat -p 6 --solexa1.3-quals -o Indx15_tophat-G -G ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf ~/iGenome/Mus_musculus/UCSC/mm9/Sequence/BowtieIndex/genome ../Indx15.txt &
nohup tophat -p 6 --solexa1.3-quals -o Indx14_tophat-G -G ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf ~/iGenome/Mus_musculus/UCSC/mm9/Sequence/BowtieIndex/genome ../Indx14.txt &
nohup tophat -p 6 --solexa1.3-quals -o Indx13_tophat-G -G ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf ~/iGenome/Mus_musculus/UCSC/mm9/Sequence/BowtieIndex/genome ../Indx13.txt &
nohup tophat -p 6 --solexa1.3-quals -o Indx12_tophat-G -G ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf ~/iGenome/Mus_musculus/UCSC/mm9/Sequence/BowtieIndex/genome ../Indx12.txt &
nohup tophat -p 6 --solexa1.3-quals -o Indx11_tophat-G -G ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf ~/iGenome/Mus_musculus/UCSC/mm9/Sequence/BowtieIndex/genome ../Indx11.txt &
nohup tophat -p 6 --solexa1.3-quals -o Indx10_tophat-G -G ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf ~/iGenome/Mus_musculus/UCSC/mm9/Sequence/BowtieIndex/genome ../Indx10.txt &
* genome-miner
Need to download UCSC dataset from iGenome to genome-miner
command: [http://genome-tech.ucsd.edu/LabNotes/index.php/Rui:RNAseq_analysis_on_HL098#Database_download]
nohup tophat -p 6 --solexa1.3-quals -o Indx3_tophat-G -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf /GenomeDB/UCSC/Mus_musculus/mm9/Sequence/BowtieIndex/genome ../Indx_seq/Indx3.txt &
Others using newly download iGenome on genome-miner (~1hr download with wget, tar 20min) start at 1:40pm
nohup tophat -p 6 --solexa1.3-quals -o Indx4_tophat-G -G ~/RNAseq/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf ~/RNAseq/iGenome/Mus_musculus/UCSC/mm9/Sequence/BowtieIndex/genome ../Indx_seq/Indx4.txt &
nohup tophat -p 6 --solexa1.3-quals -o Indx5_tophat-G -G ~/RNAseq/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf ~/RNAseq/iGenome/Mus_musculus/UCSC/mm9/Sequence/BowtieIndex/genome ../Indx_seq/Indx5.txt &
nohup tophat -p 6 --solexa1.3-quals -o Indx6_tophat-G -G ~/RNAseq/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf ~/RNAseq/iGenome/Mus_musculus/UCSC/mm9/Sequence/BowtieIndex/genome ../Indx_seq/Indx6.txt &
nohup tophat -p 6 --solexa1.3-quals -o Indx7_tophat-G -G ~/RNAseq/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf ~/RNAseq/iGenome/Mus_musculus/UCSC/mm9/Sequence/BowtieIndex/genome ../Indx_seq/Indx7.txt &
nohup tophat -p 6 --solexa1.3-quals -o Indx8_tophat-G -G ~/RNAseq/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf ~/RNAseq/iGenome/Mus_musculus/UCSC/mm9/Sequence/BowtieIndex/genome ../Indx_seq/Indx8.txt &
nohup tophat -p 6 --solexa1.3-quals -o Indx9_tophat-G -G ~/RNAseq/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf ~/RNAseq/iGenome/Mus_musculus/UCSC/mm9/Sequence/BowtieIndex/genome ../Indx_seq/Indx9.txt &
nohup tophat -p 6 --solexa1.3-quals -o Indx12_tophat-G -G ~/RNAseq/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf ~/RNAseq/iGenome/Mus_musculus/UCSC/mm9/Sequence/BowtieIndex/genome ../Indx_seq/Indx12.txt &
====cufflinks====
* Meangenemachine
nohup samtools view -h ../Indx7_tophat-G/accepted_hits.bam > accepted_hits.sam &
nohup samtools view -h ../Indx8_tophat-G/accepted_hits.bam > accepted_hits.sam &
nohup samtools view -h ../Indx9_tophat-G/accepted_hits.bam > accepted_hits.sam &
nohup samtools view -h ../Indx10_tophat-G/accepted_hits.bam > accepted_hits.sam &
nohup samtools view -h ../Indx11_tophat-G/accepted_hits.bam > accepted_hits.sam &
nohup samtools view -h ../Indx13_tophat-G/accepted_hits.bam > accepted_hits.sam &
nohup samtools view -h ../Indx14_tophat-G/accepted_hits.bam > accepted_hits.sam &
nohup samtools view -h ../Indx15_tophat-G/accepted_hits.bam > accepted_hits.sam &
cd ~/Hiseq111005/Indx7/Indx7_cufflinks
cd ~/Hiseq111005/Indx8/Indx8_cufflinks
cd ~/Hiseq111005/Indx9/Indx9_cufflinks
cd ~/Hiseq111005/Indx10/Indx10_cufflinks
cd ~/Hiseq111005/Indx11/Indx11_cufflinks
cd ~/Hiseq111005/Indx13/Indx13_cufflinks
cd ~/Hiseq111005/Indx14/Indx14_cufflinks
cd ~/Hiseq111005/Indx15/Indx15_cufflinks
nohup cufflinks -g ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf accepted_hits.sam &
* Genome-miner
mkdir ~/RNAseq//Hiseq111005/Indx3/Indx3_cufflinks
cd ~/RNAseq//Hiseq111005/Indx3/Indx3_cufflinks
nohup samtools view -h ../Indx3_tophat-G/accepted_hits.bam > accepted_hits.sam &
mkdir ~/RNAseq//Hiseq111005/Indx4/Indx4_cufflinks
cd ~/RNAseq//Hiseq111005/Indx4/Indx4_cufflinks
nohup samtools view -h ../Indx4_tophat-G/accepted_hits.bam > accepted_hits.sam &
mkdir ~/RNAseq//Hiseq111005/Indx5/Indx5_cufflinks
cd ~/RNAseq//Hiseq111005/Indx5/Indx5_cufflinks
nohup samtools view -h ../Indx5_tophat-G/accepted_hits.bam > accepted_hits.sam &
mkdir ~/RNAseq//Hiseq111005/Indx6/Indx6_cufflinks
cd ~/RNAseq//Hiseq111005/Indx6/Indx6_cufflinks
nohup samtools view -h ../Indx6_tophat-G/accepted_hits.bam > accepted_hits.sam &
mkdir ~/RNAseq//Hiseq111005/Indx12/Indx12_cufflinks
cd ~/RNAseq//Hiseq111005/Indx6/Indx6_cufflinks
nohup samtools view -h ../Indx6_tophat-G/accepted_hits.bam > accepted_hits.sam &
nohup ~/bin/cufflinks-1.1.0.Linux_x86_64/cufflinks -g ~/RNAseq/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf accepted_hits.sam &
=====comparison of genes.fpkm_tracking=====
* count cufflinks genes ID:
cd ../Indx3_cufflinks/
cd ../Indx4_cufflinks/
cd ../Indx5_cufflinks/
cd ../Indx6_cufflinks/
cd ../Indx7_cufflinks/
cd ../Indx8_cufflinks/
cd ../Indx9_cufflinks/
cd ../Indx10_cufflinks/
cd ../Indx11_cufflinks/
cd ../Indx12_cufflinks/
cd ../Indx13_cufflinks/
cd ../Indx14_cufflinks/
cd ../Indx15_cufflinks/
awk '{print $5}' genes.fpkm_tracking | sort | uniq | wc -l > count_genes.fpkm_tracking
less count_genes.fpkm_tracking
* table
There are no much difference in table organization or genes calling
The only difference in gene numbers, I believe, is due to the latest gtf I downloaded in genome-minor which leads to more annotations.
It actually means the matrix problem I met last time is not due to cufflinks files, but the script itself.
{| {{table}} border=1
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|'''tophat w/o G'''
| align="center" style="background:#f0f0f0;"|'''tophat w G'''
| align="center" style="background:#f0f0f0;"|''''''
|-
| cufflinks_Indx3||21266||21913||genome-miner
|-
| cufflinks_Indx4||21250||21896||genome-miner
|-
| cufflinks_Indx5||21238||21881||genome-miner
|-
| cufflinks_Indx6||21253||21896||genome-miner
|-
| cufflinks_Indx7||21348||21344||meangenemachine
|-
| cufflinks_Indx8||21210||21212||meangenemachine
|-
| cufflinks_Indx9||21198||21196||meangenemachine
|-
| cufflinks_Indx10||21209||21221||meangenemachine
|-
| cufflinks_Indx11||21200||21208||meangenemachine
|-
| cufflinks_Indx12||21155||21805||genome-miner
|-
| cufflinks_Indx13||21290||21292||meangenemachine
|-
| cufflinks_Indx14||21257||21255||meangenemachine
|-
| cufflinks_Indx15||21158||21159||meangenemachine
|}
====Combine_geneTracking.pl====
* Original script
[[Media:Combine_geneTracking.pl.jpg]]
* Problem 1: L27 - field[0] is the tracking_id (most of them are CUFF.xxxx) in genes.fpkm_tracking file
# Tracking_id doesn't includ all genes in "gene_short_name" field[4], for some reason.
# Same tracking_id (eg. CUFF.1) in different cufflinks files is referred to different sequence (different loci), which can't be compared!
# Change to filed[4], saved as f4_Combine_geneTracking.pl under ruiliu@genome-miner:~/RNAseq/Hiseq111005/cufflinks_UCSC
# Run new script with old cufflinks files, resulting in 7,592 genes (3,272 lines, including CUFF.xxx and genes, in old gene_expression_matrix.txt)
genes.fpkm_tracking file example:
{| {{table}}
| align="center" style="background:#f0f0f0;"|'''tracking_id'''
| align="center" style="background:#f0f0f0;"|'''class_code'''
| align="center" style="background:#f0f0f0;"|'''nearest_ref_id'''
| align="center" style="background:#f0f0f0;"|'''gene_id'''
| align="center" style="background:#f0f0f0;"|'''gene_short_name'''
| align="center" style="background:#f0f0f0;"|'''tss_id'''
| align="center" style="background:#f0f0f0;"|'''locus'''
| align="center" style="background:#f0f0f0;"|'''length'''
| align="center" style="background:#f0f0f0;"|'''coverage'''
| align="center" style="background:#f0f0f0;"|'''status'''
| align="center" style="background:#f0f0f0;"|'''FPKM'''
| align="center" style="background:#f0f0f0;"|'''FPKM_conf_lo'''
| align="center" style="background:#f0f0f0;"|'''FPKM_conf_hi'''
|-
| Xkr4||-||-||Xkr4||Xkr4||TSS1758||chr1:3204562-3661579||-||-||OK||0||0||0
|-
| CUFF.1||-||-||CUFF.1||-||-||chr1:3363176-3363727||-||-||OK||0.928728||0.449135||1.40832
|-
| CUFF.2||-||-||CUFF.2||-||-||chr1:3363898-3364505||-||-||OK||0.584168||0.231902||0.936435
|-
| CUFF.3||-||-||CUFF.3||-||-||chr1:3620867-3620976||-||-||OK||3537.97||2547.14||4528.8
|-
| CUFF.4||-||-||CUFF.4||-||-||chr1:3624844-3626291||-||-||OK||0.52528||0.336594||0.713966
|-
| CUFF.5||-||-||CUFF.5||Rp1||TSS14156||TSS17089||chr1:4280926-4399322||-||-||OK||0||0
|-
| Sox17||-||-||Sox17||Sox17||TSS20429||chr1:4481008-4486494||-||-||OK||0.0931688||0.0414881||0.14485
|-
| CUFF.6||-||-||CUFF.6||-||-||chr1:4496126-4496932||-||-||OK||0.458547||0.204191||0.712903
|-
| CUFF.7||-||-||CUFF.7||-||-||chr1:4514574-4515153||-||-||OK||1.25781||0.721481||1.79415
|-
| CUFF.8||-||-||CUFF.8||-||-||chr1:4515212-4515879||-||-||OK||5.50676||4.49926||6.51425
|}
* Problem 2: L27 - fileds[10]>0
# log(fpkm) eliminates many genes show 0 expression in early stage or certain samples
# ~half of the candidate genes showing in spermatogenesis pathway of David functional enrichment table has no match in matrix table
# Reduce to E13.5 samples, wt vs KO. 3 out of 18 don't have fpkm value in the matrix table --- is it possible to use absolute fpkm value in matrix?
# (Problem solved by Arthurva! There was a space in front of the gene name...) VLOOKUP function doesn't work between matrix table and gene list in David pathways, even after case change of all genes.
===Cufflinks w/ G annotation===
* Previously, I used -g option which allows reference as the guide
* Dr. Zhang used -G option which excludes any novel transcripts
* To repeat cufflinks in genome-miner by using the following command
nohup /home/kunzhang/softwares/cufflinks-1.1.0.Linux_x86_64/cufflinks -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf -N -p 2 -o Indx13 /home/ruiliu/RNAseq/Hiseq111005/tophat/Indx13_accepted_hits.bam > Indx13.log&
Problem need to solve:
# can .bam file directly used for cufflinks?
# is the tracking_id exactly same as gene_id?
--- genome-minor, bam file worked!
mkdir ~/RNAseq/Hiseq111005/cufflinks-G
nohup /home/kunzhang/softwares/cufflinks-1.1.0.Linux_x86_64/cufflinks -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf -N -p 2 -o Indx13 /home/ruiliu/RNAseq/Hiseq111005/tophat/Indx13_accepted_hits.bam > Indx13.log&
nohup /home/kunzhang/softwares/cufflinks-1.1.0.Linux_x86_64/cufflinks -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf -N -p 2 -o Indx14 /home/ruiliu/RNAseq/Hiseq111005/tophat/Indx14_accepted_hits.bam > Indx14.log&
nohup /home/kunzhang/softwares/cufflinks-1.1.0.Linux_x86_64/cufflinks -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf -N -p 2 -o Indx15 /home/ruiliu/RNAseq/Hiseq111005/tophat/Indx15_accepted_hits.bam > Indx15.log&
--- meangenemachine, bam file failed again!


nohup cufflinks -G ~/iGenome/ ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf -N -p 2 -o Indx13_cufflinks ~/Hiseq111005/Indx13/Indx13_tophat-G/accepted_hits.bam > Indx13.log&
==Tophat against Ensemble on 11.14.11==
*[[Rui:Tophat_11.14.11|Rui:Tophat_11.14.11]]
*[[Rui:Tophat_11.21.11|Rui:Tophat_11.21.11]]

Latest revision as of 18:07, 22 November 2011

RNAseq analysis on Hiseq111005[edit]

Comparison of program[edit]

' ' Meangenemachine Genome-miner ' '
Tophat v1.3.1 v1.3.1
Bowtie 0.12.7.0 0.12.7.0
Samtools 0.1.17 0.1.12a
BowtieIndex (UCSC) (UCSC) (Ensemble)
.fa 2011.8.8 2011.8.6 2011.10.27 2011.11.21
.ebwt 2011.5.16 2011.5.16 2011.9.27 2011.9.27
Annotation gene.gtf 2011.5.17 2011.5.17 2011.9.27 2011.9.27
file ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/gene.gtf /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf ~/RNAseq/iGenome/UCSC/Mus_musculus/UCSC/mm9/Annotation/Genes/gene.gtf [1]
cufflinks v1.0.3 v1.1.0 (2699)
cuffdiff v1.0.3 v1.1.0 (2699)

Data[edit]

  • make shortcut for input data:
ln -s /media/Ext6T/111005_SN853/RNAseq Hiseq111005/
  • make a shorcut for output data:
ln -s /media/Ext4T/DataDrive.backup/RL_Scratch/RNAseq RNAseq/
  • make copy of sequence file
cp Hiseq111005/* RNAseq/
  • combine different files in different lanes into one file for each index:
less *Indx3.txt > Indx3.txt
  • Transfer files to meangenemachine, as genome-miner is currently busy
ssh ruiliu@meangenemachine.dynamic.ucsd.edu
scp ruiliu@genome-miner:~/RNAseq/Hiseq111005/Indx_seq/* ./

Tophat_10.25.11 w/o G correction[edit]

Tophat_10.27.11 w/ G correction[edit]

Repeat RNAseq analysis on Hiseq111005 from 11.10.11[edit]


Tophat against Ensemble on 11.14.11[edit]