Rui:RNAseq analysis on Hiseq111005: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>RuiLiu
>RuiLiu
 
(38 intermediate revisions by the same user not shown)
Line 1: Line 1:
==RNAseq analysis on Hiseq111005==
==RNAseq analysis on Hiseq111005==


===Data===
===Comparison of program===
* RNA libraries were used for HL101 run, Indx3-15, [http://genome-tech.ucsd.edu/LabNotes/index.php/Rui_Liu#RNA_libraries_for_Dr._Yi_Zhang_lab
* Libraries to focus on is Indx9|10 (wt_m), Indx11|12 (wt_f), Indx13 (KO_m), Indx14|15 (KO_f)
* Data was stored in genome-miner:/media/Ext6T/111005_SN853/RNAseq
 
* 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 mapping====
=====genome-miner=====
* PATH issue from Athurva's note:
* Need PATH for bowtie, samtools, as well as GenomeDB??
 
ruiliu@genome-miner:~/RNAseq/Hiseq111005/tophat$ ~kunzhang/softwares/bowtie-latest/
bowtie                bowtie-build-debug    bowtie-inspect        doc/                  indexes/              scripts/             
bowtie-build          bowtie-debug          bowtie-inspect-debug  genomes/              reads/               
ruiliu@genome-miner:~/RNAseq/Hiseq111005/tophat$ bowtie
^C
ruiliu@genome-miner:~/RNAseq/Hiseq111005/tophat$ ls
nohup.out  tophat_Mm_Indx3
ruiliu@genome-miner:~/RNAseq/Hiseq111005/tophat$ echo $PATH
/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games
ruiliu@genome-miner:~/RNAseq/Hiseq111005/tophat$ PATH=/home/kunzhang/softwares/
abyss-1.2.5/                            Dindel/                                QuEST_2.4/
audy-stitch-db9e338/                    dindel-1.01-linux.tar.gz                QuEST_2.4.tar.gz
audy-stitch-db9e338.tar.gz              dindel-1.01-python/                    README.txt
beagle.jar                              GenomeAnalysisTK-1.0.3864/              samtools-0.1.12a/
blast-2.2.20/                          GenomeAnalysisTK-1.0.4905/              samtools-0.1.7_x86_64-linux/
blast-2.2.20-x64-linux.tar.gz          GenomeAnalysisTK-1.0.5083/              samtools-0.1.7_x86_64-linux.tar
blat_34/                                GenomeAnalysisTK-latest/                samtools-latest/
bowtie-0.12.7/                          GenomeAnalysisTK-latest.tar            SegSeq_1.0.1/
bowtie-latest/                          gm_key_64.tar                          SegSeq_1.0.1.tar.gz
brat-1.2.2/                            greatTools/                            SHERA/
brat-1.2.2.tar.gz                      greatTools.tar.gz                      SHERA_files.tar.gz
bwa-0.5.8c/                            Homo_sapiens_UCSC_hg19.tar.gz          SNVMix2-0.11.8-r3/
bwa-0.5.9/                              human_empty.bam                        SNVMix2-0.11.8-r3.tar.gz
bwa-0.5.9.tar                          impute_v2.1.2_x86_64_static.tgz        soap2.20release/
bwa-latest/                            jksrc/                                  SOAPdenovo_Release1.04/
cgatools-1.3.0.9-docs/                  jksrc.zip                              sratoolkit.2.1.6-centos_linux64/
cgatools-1.3.0.9-docs.tar.gz            macs_1.4.1.deb                          sratoolkit.2.1.6-centos_linux64.tar.gz
cgatools-1.3.0.9-linux-x86_64/          metaGene/                              stampy-1.0.8/
cgatools-1.3.0.9-linux-x86_64.tar.gz    MetaGeneMark_linux64/                  stampy-latest.tgz
cnver-0.7.2/                            MetaGeneMark_linux64.tar.gz            taoliu-MACS-7268e40/
cnver-0.7.2.tar.gz                      mga_ia64.tar                            taoliu-MACS-v2.0.7-11-g7268e40.tar.gz
CNVnator/                              ncbi-blast-2.2.24+/                    tophat-1.2.0.Linux_x86_64/
CNVnator_v0.2.2.zip                    OLB-1.8.0/                              tophat-1.3.1.Linux_x86_64/
cufflinks-1.0.3.Linux_x86_64/          OLB-1.9.3/                              tophat-1.3.1.Linux_x86_64.tar.gz
cufflinks-1.0.3.Linux_x86_64.tar.gz    OLB-1.9.3.tar.tar.gz                    tophat-latest/
cufflinks-1.1.0.Linux_x86_64/          Phrap/                                  trinityrnaseq_r2011-07-13/
cufflinks-1.1.0.Linux_x86_64.tar.gz    Phrap-distrib.tar.Z                    trinityrnaseq_r2011-07-13.tgz
cufflinks-latest/                      phred-dist-020425.c-acd.tar.Z          velvet_1.0.18/
Cython-0.15/                            picard-tools-1.38/                     
Cython-0.15.tar.gz                      picard-tools-latest/                   
ruiliu@genome-miner:~/RNAseq/Hiseq111005/tophat$ PATH=/home/kunzhang/softwares/^C
ruiliu@genome-miner:~/RNAseq/Hiseq111005/tophat$ PATH=/home/kunzhang/softwares/bowtie-latest:$PATH /home/kunzhang/softwares/tophat-latest/tophat -p 8 --solexa1.3-quals -o tophat_Mm_Indx3 -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf /GenomeDB/UCSC/Mus_musculus/mm9/Sequence/BowtieIndex/genome ../Indx_seq/Indx3.txt
 
* Add both bowtie and samtools to PATH:
PATH=/home/kunzhang/softwares/bowtie-latest:$PATH PATH=/home/kunzhang/softwares/samtools-latest:$PATH /home/kunzhang/softwares/tophat-latest/tophat -p 6 --solexa1.3-quals -o tophat_Indx9 /GenomeDB/UCSC/Mus_musculus/mm9/Sequence/BowtieIndex/genome ../Indx_seq/Indx9.txt
 
* Dr. Zhang fixed the problem of PATH, then simply run tophat under my directory:
tophat -p 6 --solexa1.3-quals -o tophat_Indx9 /GenomeDB/UCSC/Mus_musculus/mm9/Sequence/BowtieIndex/genome ../Indx_seq/Indx9.txt
 
=====meangenemachine=====
# tophat with G correction for major mRNA, based on data analysis on HL098, mapped reads from G correction are more (~100k) than ones w/o G correction
# Tophat without G correction for non-coding RNAs, plan as follows:
Here is my recommendation for your analysis:
(1)    Perform tophat mapping without any gene model.
(2)    Perform cuffdiff analysis using the UCSC gene model (like you did before) to look at protein coding genes.
(3)    Perform cuffdiff analysis using the Ensembl gene model to look at both coding and noncoding genes (you can compare the coding ones between (2) and (3) to check the consistency.
(4)    For functional annotation of LincRNAs that we don’t know too much about, perhaps you can use the GREAT analysis (http://great.stanford.edu/public/cgi-bin/greatWeb.php) , because most of the LincRNAs act in a cis- manner.
 
Ensemble reference
/GenomeDB/MmGenome/Mus_musculus.NCBIM37.64.gtf
 
mkdir tophat_Indx15
nohup tophat -p 6 --solexa1.3-quals ~/iGenome/Mus_musculus/UCSC/mm9/Sequence/BowtieIndex/genome ../Indx15.txt &
mkdir tophat_Indx15-g
nohup tophat -p 6 --solexa1.3-quals -G ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf ~/iGenome/Mus_musculus/UCSC/mm9/Sequence/BowtieIndex/genome ../Indx15.txt &
 
 
=====bowtie.left_kept_reads.fixmap.log=====
* As proxy of mapping rate.
* Two issues: 1. reads spinning multiple exons can not be included 2. clonal reads included
* Total alignment can be calculated:
samtools flagstat accepted_hits.bam
awk '{print $1}' accepted_hits.bam | wc -l
* Total reads can be estimated as
awk '{print $1}' accepted_hits.bam | sort | uniq | wc -l
* However, total reads reach to 95% of total processed reads???
 
{| {{table}} border=1
{| {{table}} border=1
| 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;"|'''bowtie.left_kept_reads.fixmap.log'''
| align="center" style="background:#f0f0f0;"|'''Meangenemachine'''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|'''Genome-miner'''
| 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;"|'''report.log'''
| align="center" style="background:#f0f0f0;"|''''''
| align="center" style="background:#f0f0f0;"|''''''
|-
|-
| ||||processed||aligned reads ||percentage||failed reads||reads sup. -m||happy splice reads||percentage
| Tophat||||v1.3.1||v1.3.1||||
|-
|-
| E13.5 wt_m1||Indx5||32,990,203||20,649,504||0.6259||12,115,555||225,144||134,070||0.6300
| Bowtie||||0.12.7.0||0.12.7.0||||
|-
|-
| E13.5 wt_f1||Indx6||27,141,595||16,807,921||0.6193||10,178,028||155,646||139,151||0.6244
| Samtools||||0.1.17||0.1.12a||||
|-
|-
| E13.5 wt_m1||Indx9||46,375,424||29,145,806||0.6285||16,865,642||363,976||143,742||0.6316
| BowtieIndex||||(UCSC)||(UCSC)||||(Ensemble)
|-
|-
| E13.5 wt_m2||Indx10||36,895,622||23,534,232||0.6379||13,046,273||315,117||144,857||0.6418
| ||.fa||2011.8.8||2011.8.6||2011.10.27||2011.11.21
|-
|-
| E13.5 wt_f1||Indx11||41,338,700||26,193,747||0.6336||14,856,837||288,116||148,068||0.6372
| ||.ebwt||2011.5.16||2011.5.16||2011.9.27||2011.9.27
|-
|-
| E13.5 wt_f2||Indx12||50,183,035||32,476,668||0.6472||17,347,409||358,958||152,887||0.6502
| Annotation||gene.gtf||2011.5.17||2011.5.17||2011.9.27||2011.9.27
|-
|-
| E13.5 KO_m1||Indx13||34,342,309||22,409,676||0.6525||11,660,995||271,638||133,327||0.6564
| ||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]
|-
|-
| E13.5 KO_f1||Indx14||40,773,127||26,933,424||0.6606||13,505,587||334,116||135,519||0.6639
| cufflinks||||v1.0.3||v1.1.0 (2699)||||
|-
|-
| E13.5 KO_f2||Indx15 ||44,644,927||29,356,424||0.6576||14,891,518||396,985||147,430||0.6609
| cuffdiff||||v1.0.3||v1.1.0 (2699)||||
|}
|}


====Transfer and convert files====
===Data===
=====Transfer files=====
* RNA libraries were used for HL101 run, Indx3-15, [http://genome-tech.ucsd.edu/LabNotes/index.php/Rui_Liu#RNA_libraries_for_Dr._Yi_Zhang_lab
* Transfer files to genome-miner from meangenemachine
* Libraries to focus on is Indx9|10 (wt_m), Indx11|12 (wt_f), Indx13 (KO_m), Indx14|15 (KO_f)
Transfer files to ruiliu@genome-miner:~/RNAseq/Hiseq111005/tophat:
* Data was stored in genome-miner:/media/Ext6T/111005_SN853/RNAseq
scp ruiliu@meangenemachine.dynamic.ucsd.edu:~/Hiseq111005/tophat_Indx15/tophat_out/accepted_hits.bam ./Indx15_accepted_hits.bam
scp ruiliu@meangenemachine.dynamic.ucsd.edu:~/Hiseq111005/tophat_Indx14/tophat_out/accepted_hits.bam ./Indx14_accepted_hits.bam
scp ruiliu@meangenemachine.dynamic.ucsd.edu:~/Hiseq111005/tophat_Indx13/tophat_out/accepted_hits.bam ./Indx13_accepted_hits.bam
scp ruiliu@meangenemachine.dynamic.ucsd.edu:~/Hiseq111005/tophat_Indx12/tophat_out/accepted_hits.bam ./Indx12_accepted_hits.bam
scp ruiliu@meangenemachine.dynamic.ucsd.edu:~/Hiseq111005/tophat_Indx11/tophat_out/accepted_hits.bam ./Indx11_accepted_hits.bam
scp ruiliu@meangenemachine.dynamic.ucsd.edu:~/Hiseq111005/tophat_Indx10/tophat_out/accepted_hits.bam ./Indx10_accepted_hits.bam
scp ruiliu@meangenemachine.dynamic.ucsd.edu:~/Hiseq111005/tophat_Indx5/tophat_out/accepted_hits.bam ./Indx5_accepted_hits.bam


* Transfer files from genome-miner to meangenemachine
* 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/* ./


ruiliu@meangenemachine-desktop:~/Hiseq111005/tophat_Indx9$
===Tophat_10.25.11 w/o G correction===
scp ruiliu@genome-miner:~/RNAseq/Hiseq111005/tophat/tophat_Indx9/* ./
*[[Rui:Tophat mapping|Tophat mapping]]
scp ruiliu@genome-miner:~/RNAseq/Hiseq111005/tophat/tophat_Indx6/* ./
*[[Rui:Transfer and convert files|Transfer and convert files]]
*[[Rui:Cuffdiff|Cuffdiff]]
*[[Rui:Cufflinks_10.26.11|Cufflinks_10.26.11]]
*[[Rui:Matrix|Matrix]]


=====Convert files=====
===Tophat_10.27.11 w/ G correction===
* Convert bam to sam in genome-miner
*[[Rui:Tophat mapping 10.27.11|Tophat mapping 10.27.11]]
*[[Rui:Cufflinks 10.28.11|Cufflinks 10.28.11]]
*[[Rui:Matrix 10.30.11|Matrix 10.30.11]]
*[[Rui:Cufflinks_11.7.11|Cufflinks_11.7.11 w/ G annotation]]


samtools view -h Indx10_accepted_hits.bam > Indx10_accepted_hits.sam
==Repeat RNAseq analysis on Hiseq111005 from 11.10.11==
samtools view -h Indx11_accepted_hits.bam > Indx11_accepted_hits.sam
* Lesson 1: Use exactly same VERSION of programs/softwares to do analysis on ALL data in ONE computer
nohup samtools view -h Indx12_accepted_hits.bam > Indx12_accepted_hits.sam &
* Lesson 2: Random check for comparison
nohup samtools view -h Indx13_accepted_hits.bam > Indx13_accepted_hits.sam
*[[Rui:Tophat_11.10.11|Tophat_11.10.11]]
nohup samtools view -h Indx14_accepted_hits.bam > Indx14_accepted_hits.sam
*[[Rui:Cufflinks_11.11.11|Cufflinks_11.11.11]]
nohup samtools view -h Indx15_accepted_hits.bam > Indx15_accepted_hits.sam
*[[Rui:Cuffdiff_11.11.11|Cuffdiff_11.11.11]]
nohup samtools view -h Indx5_accepted_hits.bam > Indx5_accepted_hits.sam
nohup samtools view -h Indx6_accepted_hits.bam > Indx6_accepted_hits.sam
nohup samtools view -h Indx9_accepted_hits.bam > Indx9_accepted_hits.sam


* Convert bam to sam in meangenemachine
nohup samtools view -h ./tophat_Indx10/tophat_out/accepted_hits.bam > ./tophat_Indx10/tophat_out/accepted_hits.sam &
nohup samtools view -h ./tophat_Indx11/tophat_out/accepted_hits.bam > ./tophat_Indx11/tophat_out/accepted_hits.sam &
nohup samtools view -h ./tophat_Indx12/tophat_out/accepted_hits.bam > ./tophat_Indx12/tophat_out/accepted_hits.sam &
nohup samtools view -h ./tophat_Indx13/tophat_out/accepted_hits.bam > ./tophat_Indx13/tophat_out/accepted_hits.sam &
nohup samtools view -h ./tophat_Indx14/tophat_out/accepted_hits.bam > ./tophat_Indx14/tophat_out/accepted_hits.sam &
nohup samtools view -h ./tophat_Indx15/tophat_out/accepted_hits.bam > ./tophat_Indx15/tophat_out/accepted_hits.sam &
nohup samtools view -h ./tophat_Indx5/tophat_out/accepted_hits.bam > ./tophat_Indx5/tophat_out/accepted_hits.sam &
nohup samtools view -h ./tophat_Indx6/tophat_out/accepted_hits.bam > ./tophat_Indx6/tophat_out/accepted_hits.sam &
nohup samtools view -h ./tophat_Indx9/tophat_out/accepted_hits.bam > ./tophat_Indx9/tophat_out/accepted_hits.sam &
====Cuffdiff====
* Cuffdiff comparison in Ensemble (genome-miner):
* /GenomeDB/MmGenome/Mus_musculus.NCBIM37.64.gtf only differs from UCSC gene.gtf in Chr. ID (eg. 8 vs Chr8)
* Dr. Zhang ran a script to change to Chr.ID
cd cuffdiff_Ensemble
cuffdiff /GenomeDB/MmGenome/Mus_musculus.NCBIM37.64.gtf -N -p 8 -o cuffDiff_wtE13.5f_vs_KOE13.5f ../topha/Indx6_accepted_hits.sam,../tophat/Indx11_accepted_hits.sam,../tophat/Indx12_accepted_hits.sam  ../tophat/Indx14_accepted_hits.sam,../tophat/Indx15_accepted_hits.sam
cuffdiff /GenomeDB/MmGenome/Mus_musculus.NCBIM37.64.gtf -N -p 8 -o cuffDiff_wtE13.5m_vs_KOE13.5m ../tophat/Indx5_accepted_hits.sam,../tophat/Indx9_accepted_hits.sam,../tophat/Indx10_accepted_hits.sam  ../tophat/Indx13_accepted_hits.sam
* Cuffdiff comparison in UCSC (meangenemachine):
cd cuffdiff_UCSC
cuffdiff ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf -N -p 8 -o cuffDiff_wtE13.5f_vs_KOE13.5f ../tophat_Indx6/tophat_out/accepted_hits.bam,../tophat_Indx11/tophat_out/accepted_hits.bam,../tophat_Indx12/tophat_out/accepted_hits.bam ../tophat_Indx14/tophat_out/accepted_hits.bam,../tophat_Indx15/tophat_out/accepted_hits.bam
cuffdiff ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf -N -p 8 -o cuffDiff_wtE13.5m_vs_KOE13.5m ../tophat_Indx5/tophat_out/accepted_hits.bam,../tophat_Indx9/tophat_out/accepted_hits.bam,../tophat_Indx10/tophat_out/accepted_hits.bam ../tophat_Indx13/tophat_out/accepted_hits.bam
* Transfer files in meangenemachine to laptop
scp ruiliu@meangenemachine.dynamic.ucsd.edu:~/Hiseq111005/cuffdiff_UCSC/cuffDiff_wtE13.5f_vs_KOE13.5f/gene_exp.diff ./Desktop/wtE13.5f_vs_KOE13.5f_gene_exp.diff_UCSC
scp ruiliu@meangenemachine.dynamic.ucsd.edu:~/Hiseq111005/cuffdiff_UCSC/cuffDiff_wtE13.5m_vs_KOE13.5m/gene_exp.diff ./Desktop/wtE13.5m_vs_KOE13.5m_gene_exp.diff_UCSC
* Transfer files in genome-miner to laptop
scp ruiliu@genome-miner:~/RNAseq/Hiseq111005/cuffdiff_Ensemble/cuffDiff_wtE13.5f_vs_KOE13.5f/gene_exp.diff ./Desktop/wtE13.5f_vs_KOE13.5f_gene_exp.diff_ensemble
scp ruiliu@genome-miner:~/RNAseq/Hiseq111005/cuffdiff_Ensemble/cuffDiff_wtE13.5m_vs_KOE13.5m/gene_exp.diff ./Desktop/wtE13.5m_vs_KOE13.5m_gene_exp.diff_ensemble
====Cufflinks====
* genes.FPKM_tracking for each Indx is used to generate a matrix
* prepare a matrix for gene set analysis
cufflinks -g ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf -o Indx5_cufflinks ../tophat_Indx5/tophat_out/accepted_hits.sam
nohup cufflinks -g ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf -o Indx6_cufflinks ../tophat_Indx6/tophat_out/accepted_hits.sam &
nohup cufflinks -g ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf -o Indx9_cufflinks ../tophat_Indx9/tophat_out/accepted_hits.sam &
nohup cufflinks -g ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf -o Indx10_cufflinks ../tophat_Indx10/tophat_out/accepted_hits.sam &
nohup cufflinks -g ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf -o Indx11_cufflinks ../tophat_Indx11/tophat_out/accepted_hits.sam &
nohup cufflinks -g ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf -o Indx12_cufflinks ../tophat_Indx12/tophat_out/accepted_hits.sam &
nohup cufflinks -g ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf -o Indx13_cufflinks ../tophat_Indx13/tophat_out/accepted_hits.sam &
nohup cufflinks -g ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf -o Indx14_cufflinks ../tophat_Indx14/tophat_out/accepted_hits.sam &
nohup cufflinks -g ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf -o Indx15_cufflinks ../tophat_Indx15/tophat_out/accepted_hits.sam &
nohup cufflinks -g ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf -o Indx3_cufflinks ../tophat_Indx3/tophat_out/accepted_hits.sam &
nohup cufflinks -g ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf -o Indx4_cufflinks ../tophat_Indx4/tophat_out/accepted_hits.sam &
nohup cufflinks -g ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf -o Indx7_cufflinks ../tophat_Indx7/tophat_out/accepted_hits.sam &
nohup cufflinks -g ~/iGenome/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf -o Indx8_cufflinks ../tophat_Indx8/tophat_out/accepted_hits.sam &
* Transfer to genome-miner
nohup scp -r ruiliu@meangenemachine.dynamic.ucsd.edu:~/Hiseq111005/cufflinks_UCSC ./cufflinks_UCSC
====Matrix====
* Modify Dr. Zhang's script by replacing folder names (such as “tophat_Mm_Indx3” -> “Indx3_cufflinks)
* Under cufflinks_UCSC folder, nano and paste modified script in, save "Combine_geneTracking.pl" (automatically save as -rw-r--r-- in current directory)
* Executive: chmod 755 Combine_geneTracking.pl
* Run: ./Combine_geneTracking.pl > gene_expression_matrix.txt
Problem: I ran tophat w/o G correction, so that output of cufflinks (genes.FPKM_tracking) has no gene annotation. Have to repeat tophat with 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: [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
==Tophat against Ensemble on 11.14.11==
# log(fpkm) eliminates many genes show 0 expression in early stage or certain samples
*[[Rui:Tophat_11.14.11|Rui:Tophat_11.14.11]]
# ~half of the candidate genes showing in spermatogenesis pathway of David functional enrichment table has no match in matrix table
*[[Rui:Tophat_11.21.11|Rui:Tophat_11.21.11]]
# 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.

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]