今天开始正式进入RNA-seq分析啦~
软件安装
由于昨天就已经创建好了新的conda的环境,所以今天只要把用到的软件都装一下就好了。
把要安装的软件分一下类型:
数据格式转换
-
sra-tools
: 这个软件用于把NCBI里下载的.sra
格式的数据给转换成fastq。用到的工具是里面的fastq-dump
数据质控
数据质量的好坏是数据分析中的关键,所谓“garbage in, garbage out
”,如果一开始的数据是有问题的,后续再怎么分析都是徒劳。
那数据会存在什么样的问题呢?
主要是二代测序增加进去adapters序列,和测序中的低质量数据。它们会影响后续的比对、定量准确性,更会影响下游分析,因此要先生成质量报告了解一下数据的质量,如果质量不好的话要把不符合要求的序列给过滤掉。
所以这一块分两步,1. 生成质量报告
;2. 过滤不符合要求的序列
。
-
fastqc
&multiqc
: 用于生成质量报告。fastqc能对每个样品生成单独的质量报告,而multiqc则是可以把多个样本的fastqc结果整合到一个文件中,这样比较方便全面的了解和对比多个数据之间的质量情况。 -
fastp
: 由陈实富大佬写的新手友好型质控软件,只要设置输入文件和输出文件就好了,推荐给各位胖友。之前写过一篇关于这个软件的简书文章:用fastp对转录组数据做QC -
trim_galore
&trimmomatic
另外两个用于做质控的软件。为了防止质控软件本身的算法原因导致分析结果的误差,条件允许的情况下建议用不同的质控软件来处理数据进行对比,取结果好的进行下游分析。
比对(mapping)
有参转录组和无参转录组分析的本质区别就在于有没有这个“参”。参是什么呢?是reference,中文世界一般叫“参考序列”。而我研究生期间做的植物基因组,就是在做这个“参”。
比对的软件有很多种,有基于基因组比对(star
、hisat2
)、基于转录组比对(bowtie
、bwa
)、不基于比对(salmon
)
我目前用bwa比较多,别的只是知道,没怎么用过。
定量
-
subread
用的是其下的featureCounts
这个小程序。这一块我不太了解,这也是我参加这次转录组学习小组的原因~
用conda安装软件
# 激活专门用于RNA-seq分析的conda环境
conda activate rnaseq
# 安装软件
conda install fastp fastqc multiqc trim-galore trimmomatic hisat2 bowtie2 subread salmon
因为我目前人在Singapore,所以用conda下载软件还挺快的。有的小朋友习惯在conda install的时候加-y参数,这样就只要等着它装好就好了,不需要再手动去确定是否安装这些软件。但是也有个问题,它在整个安装过程中不会给任何的提示和输出,就只能看到solving environment 的状态栏在“爱的魔力转圈圈”,当安装很多软件或者网络不好的时候会特别没有安全感。
数据格式转换
我想想我一般是怎么写这个的:
#!/usr/bin/env bash
list="SRR1039509.sra
SRR1039512.sra
SRR1039513.sra
SRR1039516.sra
SRR1039517.sra
SRR1039520.sra
SRR1039521.sra"
for i in $list
do
fastq-dump --gzip --split-3 ${i} -O ../01raw
done
我的list
是通过下面的命令生成的:
ls *.sra | xargs -l
而豆豆提供的脚本是这个画风:
raw=~/rnaseq/raw # 这个路径根据个人需求进行修改即可
cat $raw/srr.ids | while read i ;do
echo $i
# time fastq-dump --gzip --split-3 -A $i $raw/${i}.sra -O $raw;
done
确实看起来比我高明得多呢。学习了学习了。用time
这个命令计算运行时间是一个很好的习惯。
先用井号注释掉了真正要运行的命令,这在写脚本的过程中是个好习惯,因为有时自己会搞错变量(比如这里的i) ,于是用echo $i可以先看看是不是真正要用的,先运行一下看看结果
把脚本写进一个叫做fastqdump.sh
的文件,然后nohup运行。
(不要随便起个名字噢,不然过段时间就忘了这个当时写了是做什么用的了呢,后患无穷,整理起来特别头疼。)
我用nohup主要是它可以保留软件运行过程中的日志和报错,当然不用nohup也可以做到这点:
bash fastqdump.sh 1>fastqdump.log 2>&1
# 或者
bash fastqdump.sh 1>fastqdump.log 2>fastqdump.err
1>fastqdump.log 表示将结果的正确日志文件输出到fastqdump.log中,2>&1 这个2表示错误日志,将它也合并到1的正确日志中。
当然,如果不想要错误信息,可以把它丢到linux"黑洞"中去: 2>/dev/null
下载参考序列
根据文章介绍,所用的reference是hg19的版本。因此就下载这个版本的genome和gtf数据备用。
wget -c http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz
wget -c ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz
太晚了。明天继续写。
Hans Chen
2019年6月9日于 Westwood Hostel, Jurong West Ave. 5, Singapore.