Rui:RNAseq analysis on Hiseq111005
Jump to navigation
Jump to search
RNAseq analysis on Hiseq111005
Data
- 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 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&