我的第一次RNAseq分析

来源：互联网发布：淘宝会员管理系统编辑：程序博客网时间：2024/05/22 20:44

从师姐那儿拿到了整个大概的protocol

0. 预处理。将公司给的fastq格式的文件变成clean.fastq。去掉一些不合格的reads

1.比对。运用以下命令生成一个参考基因组。我有师姐之前留下的基因组文件，名叫GCF_000002035.5_GRCz10_genomic.fna，我也不知道这是个什么格式。从ensemble上下载的目标基因组是.fa格式的。

cd /public/home/iemb315/Danio_genome
/public/opt/sc/program/STAR-STAR_2.4.1c/source/STAR \
--runThreadN 12 --runMode genomeGenerate --genomeDir ./ \--genomeFastaFiles ./GCF_000002035.5_GRCz10_genomic.fna

得到：

[iemb315@node100 ~/Danio_genome]$ls -shltotal 15G   0 -rw-rw-r-- 1 iemb315 iemb315    0 Mar  6 22:23 Aligned.out.sam8.0K -rw-rw-r-- 1 iemb315 iemb315 6.2K Mar  8 09:56 chrLength.txt 24K -rw-rw-r-- 1 iemb315 iemb315  22K Mar  8 09:56 chrNameLength.txt 16K -rw-rw-r-- 1 iemb315 iemb315  16K Mar  8 09:56 chrName.txt 12K -rw-rw-r-- 1 iemb315 iemb315  12K Mar  8 09:56 chrStart.txt1.3G -rw-rw-r-- 1 iemb315 iemb315 1.3G Mar  6 21:51 GCF_000002035.5_GRCz10_genomic.fna1.6G -rw-rw-r-- 1 iemb315 iemb315 1.6G Mar  8 10:30 Genome4.0K -rw-rw-r-- 1 iemb315 iemb315  538 Mar  8 09:55 genomeParameters.txt   0 -rw-rw-r-- 1 iemb315 iemb315    0 Mar  6 22:23 indexForCFpe100.log104K -rw-rw-r-- 1 iemb315 iemb315 104K Mar  8 10:30 Log.out   0 -rw-rw-r-- 1 iemb315 iemb315    0 Mar  6 22:23 Log.progress.out4.0K -rw------- 1 iemb315 iemb315  234 Mar  6 22:23 nohup.out 11G -rw-rw-r-- 1 iemb315 iemb315  11G Mar  8 10:30 SA1.5G -rw-rw-r-- 1 iemb315 iemb315 1.5G Mar  8 10:30 SAindex
然后我把GCF_000002035.5_GRCz10_genomic.fna 移动走了。事实证明它在不在也没有关系了。

2.比对。用STAR（https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4631051/）这个经典的比对工具将reads比对到目标基因组上，这个目标基因组要从ensemble上下载。我有师姐之前留下的基因组文件，名叫GCF_000002035.5_GRCz10_genomic.fna

然后就是按照protocol里的basic protocol比对了。注意：首先要进入目标基因组文件夹。我给这个突变体起名叫做march

cd /public/home/iemb315/xhb/march_genome/public/opt/sc/program/STAR-STAR_2.4.1c/source/STAR \--runThreadN 12 \--genomeDir /public/home/iemb315/Danio_genome \ #这个目录不要动，里面有莫名其妙的基因组文件，可以用来比对--readFilesIn /public/home/iemb315/xhb/M1_H7NNKALXX_L2_1.clean.fq.gz /public/home/iemb315/xhb/M1_H7NNKALXX_L2_2.clean.fq.gz \--readFilesCommand zcat \--twopassMode Basic

大概跑了一个小时。野生型也是按照这个跑的。得到2个sam格式的文件。

3. sam转换成bam文件，会变小一点点。注意samtools和STAR都是不用安装的。

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

然后还是用samtools把reads排序。这个时候突然说需要一个pysam？

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

4. 配置python运行环境。

这是最蛋疼的地方了。因为接下来要定量分析。比较RPKM值。

RPKM,Reads Per Kilobase of exon model per Million mapped reads [Mortazavi etal., 2008]。还有个词叫FPKM（fragments Per Kilobase of exon model per Million mapped reads），差不多同义。

需要使用cufflink或者HTSeq这两个python包。

然鹅，服务器上找不到这两个工具，我也没有服务器的root权限。首先我在自己的目录下安装python3，然后安装easy_install的时候总是出错。

我搞不懂怎么回事，只能等第二天问师兄。

师兄他说没安装权限的话找管理员。那个管理员老师整天闲云野鹤似的我上哪儿找去

好气啊

不过我又研究了一上午，弄好了，主要参考了这个帖子http://blog.csdn.net/dream_angel_z/article/details/51338546

其实还是很简单的。

先把python安装到自己建的文件夹下，注意一定要是python2.7，更改安装目录，不然安装到root下面了

./configure -prefix='/public/home/iemb315/python27

然后下载各种二进制的安装包，解压到python27/Python-2.7.10文件夹下。

然后安装pip，用pip 安装numpy，然后执行这个就可以了。

/public/home/iemb315/python27/Python-2.7.10/bin/python setup.py install

我不知道在我的根目录下新建一个software文件夹，能不能把这些软件安装在那儿。

我可以在我的/bin下找到下一步的安装包。我安装了htseq

/public/home/iemb315/python27/Python-2.7.10/bin/htseq-count

#总结#
配置python环境.参考：http://blog.csdn.net/dream_angel_z/article/details/51338546 #在指定目录下安装python./configure -prefix='/public/home/iemb315/python3/Python-2.7.10'makemake install
#安装setuptools,不然没法装pipcd /public/home/iemb315/python27/Python-2.7.10/ #据说这些软件可以另外建一个software文件夹？wget --no-check-certificate http://pypi.python.org/packages/source/s/setuptools/setuptools-2.0.tar.gztar -xzvf setuptools-2.0.tar.gzcd setuptools-2.0/public/home/iemb315/python27/Python-2.7.10/bin/python setup.py install
#tar的用法#-f 使用档案文件或设备，这个选项通常是必选的。#-x 从档案文件中释放文件。#-v 详细报告tar处理的文件信息。如无此选项，tar不报告文件信息。#-z 用gzip来压缩/解压缩文件，加上该选项后可以将档案文件进行压缩，但还原时也一定要使用该选项进行解压缩。#安装pipcd /public/home/iemb315/python27/Python-2.7.10/ #据说这些软件可以另外建一个software文件夹？wget --no-check-certificate https://pypi.python.org/packages/41/27/9a8d24e1b55bd8c85e4d022da2922cb206f183e2d18fee4e320c9547e751/pip-8.1.1.tar.gz#md5=6b86f11841e89c8241d689956ba99ed7tar -xzf pip-8.1.1.tar.gzcd pip-8.1.1/public/home/iemb315/python27/Python-2.7.10/bin/python setup.py install#安装numpy，matplotlib，不然HTSeq没法用cd /public/home/iemb315/python27/Python-2.7.10/bin/./pip install numpy./pip install matplotlib#用这个命令安装解压后的其他二进制软件，比如这里还安装了叫做pysam的东西/public/home/iemb315/python27/Python-2.7.10/bin/python setup.py install

5. 定量

使用python包htseq。

/public/home/iemb315/python27/Python-2.7.10/bin/htseq-count -f bam -r pos -s no -i Parent --mode=intersection-nonempty \/public/home/iemb315/xhb/wt_genome/WT_march.sorted.bam.bam \/public/home/iemb315/Danio_genome/Danio_rerio.GRCz10.87.chr.gff3 >WT.count/public/home/iemb315/python27/Python-2.7.10/bin/htseq-count -f bam -r pos -s no -i Parent --mode=intersection-nonempty \/public/home/iemb315/xhb/march_genome/MT_march.sorted.bam \/public/home/iemb315/Danio_genome/Danio_rerio.GRCz10.87.chr.gff3 >MT.count

这里跑了各1.5小时，出现服务器拥堵的问题。