Rui:RNAseq analysis on Hiseq111005

From ZhangLabWiki
Revision as of 19:01, 10 November 2011 by >RuiLiu (→‎Tophat w/o G correction)
Jump to navigation Jump to search

RNAseq analysis on Hiseq111005

Data

  • 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 w/o G correction

Tophat w/ G correction

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: [1]
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.
' tophat w/o G tophat w G '
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: 
tracking_id class_code nearest_ref_id gene_id gene_short_name tss_id locus length coverage status FPKM FPKM_conf_lo 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&