一个RNA-seq的反思 | 生信菜鸟团RNA-seq 数据的简单分析 - OA_maque - 博客园
主要介绍如何分析RNA-seq 数据
以下的文档就按上面的这个教程来组织
需要安装的软件：
sra-tools，samtools, bam-readcount, bowtie, bowtie2, tophat, star, cufflinks, htseq-count, R, cummeRbund, fastqc, picard-tools, and samstat
其中已经安装的：
以下以一个例子来说明如何进行一般的 rna-seq的分析
GEO number : GSE66666
从中了解实验是如何设计的，想解决什么问题，多少样本，该研究所发表的文章等相关信息。
下载原始序列
原始序列一般以 SRA 的格式保存在 NCBI 上。
推荐写一个 sh 脚本，批量下载，即列出所有的 链接。
然后使用 sra-tools 把 sra 格式序列转换成 fq 格式序列
脚本如下：
#!/bin/bash
cd /home/user/download/myfile/
wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP%2FSRP055%2FSRP055992/SRR1871481/SRR1871481.sra
wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP%2FSRP055%2FSRP055992/SRR1871482/SRR1871482.sra
wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP%2FSRP055%2FSRP055992/SRR1871483/SRR1871483.sra
wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP%2FSRP055%2FSRP055992/SRR1871484/SRR1871484.sra
wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP%2FSRP055%2FSRP055992/SRR1871485/SRR1871485.sra
wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP%2FSRP055%2FSRP055992/SRR1871486/SRR1871486.sra
# use sra-tools to transform
export PATH=$PATH:/home/user/sratoolkit.2.5.2-ubuntu64/bin
fastq-dump *.sra
这样就把原始序列 fq 文件得到了。
为了后面分析方便，把相应的序列文件名改成相应的组
mv SRR1871481.fastq WT_Rep1.fastq
mv SRR1871482.fastq WT_Rep2.fastq
mv SRR1871483.fastq WT_Rep3.fastq
mv SRR1871484.fastq athb1_Rep1.fastq
mv SRR1871485.fastq athb1_Rep2.fastq
mv SRR1871486.fastq athb1_Rep3.fastq
Pre-Alignment QC
使用fastqc 软件来对原始序列fastq 文件进行质量检测
export PATH=/home/maque/FastQC/:$PATH
fastqc *.fastq
这样每个 fastq 文件都能生成一个 html 报告文件，很详细
使用 tophat 和 bowtie 进行比对
在开始之前，需要下载拟南芥的基因组序列，注释文件以及 INDEX文件
选择 Ensembl tigr10 版本, 解压
cd /media/文档/rna-seq-arabi/
#原始序列与 tigr10 文件夹都放在该文件夹下
export PATH=/home/maque/samtools-1.2/bin:$PATH
export PATH=/home/maque/tophat-2.1.0.Linux_x86_64/:$PATH
export PATH=$PATH:/home/maque/bowtie-1.1.2
tophat2 -p 8 --bowtie1 -G Arabidopsis_thaliana_Ensembl_TAIR10/Arabidopsis_thaliana/Ensembl/TAIR10/Annotation/Genes/genes.gtf -o WT_Rep1_output
Arabidopsis_thaliana_Ensembl_TAIR10/Arabidopsis_thaliana/Ensembl/TAIR10/Sequence/BowtieIndex/genome WT_Rep1.fastq
# 其他5个 fastq 文件与上面一致
# --bowtie1
使用bowtie1 , 默认是bowtie2
后面跟注释文件 gtf
后跟输出文件夹
# 最后面跟 原始序列 fastq 文件
这些过程完成后，说明已经完成比对，在每个新建的文件夹里面，应该有一些信息，最主要的是生成了一个 accepted_hits.bam 文件， 这个就是 samtools 生成的，后面主要也是利用这个文件继续分析。
建议提前利用 IGV 软件查看一下 该 bam 文件，可以知道mapping 的情况。
如果想查看，需要先对 该bam文件进行 index ,比如：
samtools index
WT_Rep1_output/accepted_hits.bam
Use Cufflinks to generate expression estimates from the SAM/BAM files
export PATH=/home/maque/cufflinks-2.2.1.Linux_x86_64/:$PATH
cufflinks -p 8 -o WT_Rep1_cuffout WT_Rep1_output/accepted_hits.bam
# 其他5个与之类似
# -o WT_Rep1_cuffout
最后面跟相应的 bam 文件
该过程完成后，会生成相应的文件夹，里面有相应的文件，后面会使用 transcripts.gtf 文件
Differential Expression
ls -1 *cuffout/transcripts.gtf & assembly_GTF_list.txt
cuffmerge -p 8 -o merged -g Arabidopsis_thaliana_Ensembl_TAIR10/Arabidopsis_thaliana/Ensembl/TAIR10/Annotation/Genes/genes.gtf -s Arabidopsis_thaliana_Ensembl_TAIR10/Arabidopsis_thaliana/Ensembl/TAIR10/Sequence/WholeGenomeFasta/genome.fa assembly_GTF_list.txt
# -o merged
后跟参考基因的gtf文件
后跟基因组序列文件
最后跟上一步新建的 assembly_GTF_list.txt
接下来使用 cuffdiff
cuffdiff -o rna_de/diff_out -p 8 -L WT,athb merged/merged.gtf WT_Rep1_output/accepted_hits.bam,WT_Rep2_output/accepted_hits.bam,WT_Rep3_output/accepted_hits.bam athb_Rep1_output/accepted_hits.bam,athb_Rep2_output/accepted_hits.bam,athb_Rep3_output/accepted_hits.bam
后跟输出文件目录
# -p 8 使用8线程
# -L WT,athb
'-L' tells cuffdiff the labels to use for samples
# 后跟 上一步由 cuffmerge 生成的 merged.gtf 文件
# 最后跟6个bam 文件， 组内由逗号隔开，组间由空格隔开。
该过程会新建一个diff_out 文件夹，里面包含了很多信息
这些信息可以使用 R 包 cummeRbund 很方便的进行分析
使用cummeRbund
可以按照推荐流程文档中的步骤去分析
下面主要说一些注意点：
source(&http://bioconductor.org/biocLite.R&)
biocLite(&cummeRbund&)
首先先 cd 到上一个步骤生成的 diff_out 文件夹
library(cummeRbund)
cuff &- readCufflinks()
这样即完成读入数据。
各种分析图
可以按照推荐流程中的去分析，也可以参考
寻找差异表达基因
大部分进行RNA-Seq 实验的目的主要还是寻找实验组与对照组之间的差异表达基因。
一种方法是:
mySigGeneIds&-getSig(cuff,alpha=0.05,level='genes')
mySigGeneIds
length(mySigGeneIds)
mySigGenes&-getGenes(cuff,mySigGeneIds)
mySigGenes
diffData(mySigGenes)
featureNames(mySigGenes)
另一种方法是：
gene_diff_data &- diffData(genes(cuff))
sig_gene_data &- subset(gene_diff_data, (significant == 'yes'))
sig_gene_data
这些方法列出的 gene_id 是像 XLOC_000268 这样的格式， 它所对应的通用的gene_id 比如AT1G06080 , 这种一一对应关系文件可以在合并的 merged.gtf 文件中寻找
而AT1G06080 这种gene_id 所对应的基因类型，基因名称等信息可以在 tair10 文件夹中的 gene.gtf 文件中找到
AT1G06080 这种gene_id 所对应的基因名称也可以在 同一文件夹中的 refFlat.txt 文件中找到。
也可以先把上一步输出的 gene_id 放到一个 list.txt 中，注意不要有空行，最好使用 vim , 然后利用 grep 即可：
merged.gtf
| less - S
以上就是rna-seq 数据分析的简单过程，很多细节没有提，而且还有很多其他步骤可以进行扩展，这些还需要再进一步深入理解。更多最新文章相关作者文章搜狗：感谢您阅读RNA-seq结果怎么才能看懂?答案全在这些图里---（2）基础分析结果篇 本文版权归原作者所有，本文由网友投递产生，如有侵权请联系 ，会第一时间为您处理删除。关注微口网微博
微信号：iweikou
熟悉我的人都知道RNA-seq是我的拿手好戏（如果你不熟悉我，今天过后请记住）。
但是我今天处理了一个公共数据，比对率低的惊人。
究竟为什么会发生这种小概率事情呢？
是测序数据质量不好？
对了，顺便说一句，我已经把下好的sra数据删除了!（你想笑就笑吧）
我痛恨我自己直接拷贝了博客的脚本！我痛恨&--这样的参数设置！
因为太熟练了，所以我先想了各种复杂的情况而没有一开始就从最简单的问题入手。
这篇文章告诉你，老司机如我，也有翻车的时候。数据比较新，所以理所当然的认为测序数据肯定是OK的。
0.62% overall alignment rate
此处请脑补我内心的呐喊！
因为太熟练了，所以我把wget跟fastq-dump一起运行，很久没关心过log文件。
因为太熟练了，所以我已经很久没有关心过完整的fastqc报告。
不过，暂时看到这个问题，我就试着解决一下吧，先从这个思路来，而且比对工具里面本来就有这个选项，没必要自己来trim的！
数据下载地址
~/biosoft/HISAT/current/hisat2 -p 5 -x $reference -U SRR3589962.fastq -S Akap952.sam 2&Akap952.log
(18)一个mRNA-seq实战-生信菜鸟团博客2周年精选文章集
回过头看了看fastqc的报告，发现前面10个碱基的确有问题的！如果只是对RNA-seq进行定量，可能需要trim掉，但是我以前从来不trim，照样不影响比对。（如果你认真研究过测序流程，你应该知道前面这十个左右的碱基其实是bias）
是--不—是全角半角害死人，而且这个参数不识别它居然不停止，而是忽略我的参数！
下载sra数据，转换为fastq的过程我就不讲解了！(请翻看我以前的文章)
果然，质量值好到爆！
0.49% overall alignment rate起初我怀疑是参考基因组用错了，但是我查看了GEO里面的介绍，的确是mouse的ESC，所以我用grcm38没有问题！然后我怀疑是测序数据质量的问题，但是质量再差也不会导致如此低的比对率啊！所以我还是用fastqc检查了一下：
(20)一个WES实战-生信菜鸟团博客2周年精选文章集
终于，我想起了最开始的fastqc，决定再回过头看看测序数据的fastqc报告，请看下图，这么重要的图我居然忽略掉了，忽略掉了，略掉了，掉了，了。。。
是--不—是全角半角害死人，而且这个参数不识别它居然不停止，而是忽略我的参数！
https://ccb.jhu.edu/software/hisat2/index.shtml
唯一值得你看的就是上面这个图是--不—是全角半角害死人，而且这个参数不识别它居然不停止，而是忽略我的参数！
更要命的是我把wget跟fastq-dump一起运行，而wget会给出一大堆的log日志把fastq-dump的报错日志给掩盖了。
下面是我修改后的代码！cut -f 3 config.txt | do wget $id 2&/dev/done
貌似的确是mouse基因组的染色体长度呀！很诡异，而且我清楚的记得，我下载的就是mouse的基因的索引。这里也没有问题。
那么就是我hisat2这个步骤的问题，我首先怀疑是不是我下载hisat的index搞错了，虽然看起来我命名是grcm38，但是有可能是我下载错误！我打开了sam文件看了看开头：
(17)一个ChIP-seq实战-生信菜鸟团博客2周年精选文章集
-3/--trim3 &int& trim &int& bases from 3'/right end of reads (0)所以我加上了-p 6 -5 10 -3 10 --local 参数，比对人，可以拿到 35.60% overall alignment rate，比对mouse，可以拿到 98.80% overall alignment rate ，我勒个去，问题出来了，看起来好像是应该trim掉呀。以前的万能默认参数不行了？
spots for SRR3589962.sra我用的是hisat2工具来比对，一般情况下就用默认参数。reference=/home/jianmingzeng/reference/index/hisat/grcm38/genome
难道grcm38与mm10有差别？ 我就先用bowtie2测试一下mm10吧，毕竟我还没有hisat2的mm10的index。head -1000 SRR3589959.fastq &tmp.fq
~/biosoft/bowtie/bowtie2-2.2.9/bowtie2 -p 6 -x ~/reference/index/bowtie/mm10 -U tmp.fq -S tmp.sam结果我挑出来的这1000条序列，全军覆没了，0.00% overall alignment rate，我傻眼了！ 没办法呀，逼着我换hg19参考基因组看看：~/biosoft/bowtie/bowtie2-2.2.9/bowtie2 -p 6 -x ~/reference/index/bowtie/hg19 -U tmp.fq -S tmp.hg19.sam仍然是全军覆没了，0.00% overall alignment rate，心好累，继续傻眼！
我然后用同样的参数，测试了hisat2工具，但是hisat2里面压根就没有local的选项，仅仅是trim一下，对比对的改善毫无意义，所以重点在于--local这个参数，但它只是表象，本质还是这个测序数据出问题了！
~/biosoft/HISAT/current/hisat2 -p 5 -x $reference -U SRR3589959.fastq -S control1.sam 2&control1.log
~/biosoft/HISAT/current/hisat2 -p 5 -x $reference -U SRR3589960.fastq -S control2.sam 2&control2.log
数据为什么会出问题呢?
因为太熟练了，所以我对用了一千遍fastq-dump的命令产生了一万个放心。
然后我检查我的脚本，我自己从我的博客里面复制了我的代码。
spots for SRR3589960.sra
ls *sra | do ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --gzip --split-3 $done这就是老司机翻车以及我debug的心路历程，希望你们引以为戒！
而且我再去GEO介绍上面看，上面赫然写着PAIRED！ 我死也想不明白，我明明是加了--split-3 参数呀，为什么sra转换成fastq会出这么明显的错误呢？
因为前面一直处理的是单端的数据，所以这个错误没有被发现。
ls sam |do (nohup samtools sort -n -@ 5 -o ${id%%.}.Nsort.bam $id &);done但是让我意外的是比对率出奇的低！请看我结果。。。0.48% overall alignment rate
~/biosoft/HISAT/current/hisat2 -p 5 -x $reference -U SRR3589961.fastq -S Akap951.sam 2&Akap951.log
spots for SRR3589961.sra
但是有个问题，虽然我用local模式都比对上了，但是首先100bp的reads我切成了80，而且都是40M，40S，说明只有reads的一半成功比对到了参考基因组序列！
再联想到前面的40M，40S我瞬间明白了，这肯定是一个双端测序被我搞成了单端测序数据！
0.48% overall alignment rate
具体参数网站：https://ccb.jhu.edu/software/hisat2/manual.shtml-5/--trim5 &int& trim &int& bases from 5'/left end of reads (0)
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE81916
比对工具的默认参数不行？
所以，一切看起来都是因为我是一个老司机！
而且我抽取了几条序列去blat一下，发现也可以比对，而且明显是跨越intron的比对，超级经典的RNA-seq数据呀，这完全没问题啊。( 其实我这个blat结果也没有看仔细，正常的reads不应该被截成比对到基因组的正负链的，这其实预示着我把PE序列拼接了。)
总之，这次处理了GEO里面的4个公共数据，数据量如下：Written
spots for SRR3589959.sra
（图一十一）
难道grcm38与mm10有差别？
(19)一个affymetrix表达芯片实战-生信菜鸟团博客2周年精选文章集
如果你认为自己还是一个菜鸟就请记住老司机的这句话：行走江湖，靠的是一个稳字文：Jimmy
编辑校对：一只思考问题的熊
（图一十二）
看过本文的人还看过
人气：1334 更新：
人气：656 更新：
人气：536 更新：
人气：481 更新：
生信菜鸟团的更多文章
大家在看 ^+^
推荐阅读 ^o^
『中國邊疆研究與歷史書寫』研討會日程安排
过真伤己、过直伤人
中国人走得太远太快，灵魂跟不上了（深度好文）
他说第二，有人敢说第一吗？
猜你喜欢 ^_^
24小时热门文章
微信扫一扫
分享到朋友圈}