我的第二次RNAseq分析

数据还是上次那些，因为我没法解决.count文件里面的数字全是0的情况。所以打算从头开始做。

这次我对流程更加熟练了，另外又看了不少攻略，理解也更深了一些

##########################

1 基因组索引的建立

##########################

需要从ensemble的FTP下载一个.fa文件和一个.gtf文件（或者.gff3）

genomeFastaFiles 这个参数可以使用  .fna     .fa这两种文件格式--sjdbGTFfile这个选项可以使用.gff3   .gtf这两种格式可以用。

这次我用的是Danio_rerio.GRCz10.dna.toplevel.fa 和Danio_rerio.GRCz10.87.gtf 这两个文件。

不要用GCF_000002035.5_GRCz10_genomic.fna 和Danio_rerio.GRCz10.87.gtf的组合！

不要用Danio_rerio.GRCz10.dna.toplevel.sm.fa和Danio_rerio.GRCz10.87.gtf的组合！

一定要用这个文件啊 Danio_rerio.GRCz10.dna.toplevel.fa

命令如下：

/public/opt/sc/program/STAR-STAR_2.4.1c/source/STAR \--runThreadN 20 \--runMode genomeGenerate \--genomeDir /public/home/iemb315/genome \--genomeFastaFiles /public/home/iemb315/Danio_genome/Danio_rerio.GRCz10.dna.toplevel.fa \--sjdbGTFfile /public/home/iemb315/Danio_genome/Danio_rerio.GRCz10.87.gtf \--sjdbOverhang 149 

得到如下：

Mar 13 11:28:11 ..... Started STAR run
Mar 13 11:28:11 ... Starting to generate Genome files
Mar 13 11:29:02 ... starting to sort  Suffix Array. This may take a long time...
Mar 13 11:29:20 ... sorting Suffix Array chunks and saving them to disk...
Mar 13 11:39:33 ... loading chunks from disk, packing SA...
Mar 13 11:40:41 ... Finished generating suffix array
Mar 13 11:40:41 ... starting to generate Suffix Array index...
Mar 13 11:52:15 ..... Processing annotations GTF
Mar 13 11:52:30 ..... Inserting junctions into the genome indices
Mar 13 11:58:43 ... writing Genome to disk ...
Mar 13 11:58:45 ... writing Suffix Array to disk ...
Mar 13 11:58:59 ... writing SAindex to disk
Mar 13 11:59:01 ..... Finished successfully

如果需要根据不同基因组调整参数需要自己阅读STAR的manual

或许会总是报错 Segmentation fault (core dumped) 。管理员说是段错误。多试几次可能就好了。

##########################

2 mapping reads

##########################

把公司给的文件用filezillar 上传然后解压

gunzip /public/home/iemb315/xhb/W1_H7NNKALXX_L2_1.clean.fq.gz
gunzip /public/home/iemb315/xhb/W1_H7NNKALXX_L2_2.clean.fq.gz
gunzip /public/home/iemb315/xhb/M1_H7NNKALXX_L2_1.clean.fq.gz
gunzip /public/home/iemb315/xhb/W1_H7NNKALXX_L2_2.clean.fq.gz

解压后会覆盖掉原有的压缩文件，运行head W1_H7NNKALXX_L2_2.clean.fq 查看一条read多长

比对WT

cd /public/home/iemb315/xhb/wt/public/opt/sc/program/STAR-STAR_2.4.1c/source/STAR \--runThreadN 20 --genomeDir ~/genome \--readFilesIn ~/xhb/W1_H7NNKALXX_L2_1.clean.fq ~/xhb/W1_H7NNKALXX_L2_2.clean.fq \--sjdbGTFfile /public/home/iemb315/Danio_genome/Danio_rerio.GRCz10.87.gtf \--sjdbOverhang 149 

比对MT

cd /public/home/iemb315/xhb/mt/public/opt/sc/program/STAR-STAR_2.4.1c/source/STAR \--runThreadN 12 --genomeDir ~/genome \--readFilesIn ~/xhb/M1_H7NNKALXX_L2_1.clean.fq ~/xhb/M1_H7NNKALXX_L2_2.clean.fq \--sjdbGTFfile /public/home/iemb315/Danio_genome/Danio_rerio.GRCz10.87.gtf \--sjdbOverhang 149

在各自的目录下得到两个很大的.sam文件，以及一些运行日志。

/public/home/iemb315/xhb/wt

/public/home/iemb315/xhb/mt

###################################################

3 使用samtools预处理sam文件，这个工具不用安装，服务器本身就有。

###################################################

进入wt mt的目录

使用view

samtools view -bS Aligned.out.sam > MT_march.bam
samtools view -bS Aligned.out.sam > WT_march.bam

使用sort

samtools sort -m 10G -@ 2 WT_march.bam WT_march.sorted 
samtools sort -m 10G -@ 2 MT_march.bam MT_march.sorted

一定要查看一下生成的bam文件

samtools view -H WT_march.sorted.bam > header

samtools view -H MT_march.sorted.bam > header

#############################
4 使用htseq-count给reads计数
#############################

~/python27/Python-2.7.10/bin/htseq-count -f bam -r pos -s no --mode=intersection-nonempty ./WT_march.sorted.bam ~/Danio_genome/Danio_rerio.GRCz10.87.gtf >WT.count 

~/python27/Python-2.7.10/bin/htseq-count -f bam -r pos -s no --mode=intersection-nonempty ./MT_march.sorted.bam ~/Danio_genome/Danio_rerio.GRCz10.87.gtf >MT.count

如果所计算出的count全部为0,有可能是因为比对用的gff文件和htseq-count用的gff文件不一致造成的 现在来看是.fa文件版本不对

######################################################
5 使用DESeq2处理.count文件

######################################################

把两个.count文件下载到自己windows的环境变量目录下。

安装anaconda 安装Rstudio

DESeq2的安装

source('http://bioconductor.org/biocLite.R')
biocLite('DESeq2')

library('DESeq2')  directory<-'~/test' sampleFiles<-grep('t',list.files(directory),value=TRUE)  # sampleFiles  sampleCondition<-c('mt','wt')  sampleTable<-data.frame(sampleName=sampleFiles, fileName=sampleFiles, condition=sampleCondition) ddsHTSeq<-DESeqDataSetFromHTSeqCount(sampleTable=sampleTable, directory=directory, design=~condition)  ddsHTSeq  #######确定顺序  colData(ddsHTSeq)$condition<-factor(colData(ddsHTSeq)$condition, levels=c('mt','wt'))     dds<-DESeq(ddsHTSeq)  res<-results(dds)  res<-res[order(res$padj),]  head(res)  write.table(res, file="ddseq2.xls", sep="\t",quote=F)

                                             
        0        0           
	
				
	
					   
我的第二次RNAseq分析
	  	   
我的第一次RNAseq分析
	  	   
我的第二次建模
	  	   
我的第二次初恋
	  	   
我的第二次实验
	  	   
第二次 尝试 我的
	  	   
我的第二次面试经历
	  	   
我的C++第二次作业
	  	   
我的第二次实验命令：
	  	   
下载我完成的tinyMCE(第二次修改)
	  	   
我第二次用SecureCRT时的糗事
	  	   
我的第二次找工作之旅
	  	   
我的c++第二次程序1
	  	   
我的C++第二次实验作业
	  	   
我的C++第二次实验报告
	  	   
我的第二次上机作业项目四
	  	   
我的第二次实验，第四个项目
	  	   
第二次启动android app的过程分析
	     		  
	  	   
[BZOJ3223][Tyvj1729]文艺平衡树
	  	   
四元傅里叶显著性图-四元数-Matlab编程
	  	   
Android多线程多进程学习网址
	  	   
图像处理之简化色彩(含OpenCV代码)
	  	   
int64
	  	   
我的第二次RNAseq分析
	  	   
LCD驱动程序(四)测试
	  	   
求n以内的最大素数，若n最大为21亿
	  	   
ConcurrentHashMap复合操作问题
	  	   
做正确的事
	  	   
MySQL数据库事务隔离级别
	  	   
04. Oracle 11g 数据库关闭与启动
	  	   
HDU 2191 悼念512汶川大地震遇难同胞——珍惜现在，感恩生活【多重背包】
	  	   
触摸屏驱动之概念介绍