转录组数据分析01-比对

3 min read RNA-seq

最近有一批转录组数据需要分析，从公司拿到raw data之后进行后续分析。目前转录组测序十分成熟，测序质量都很高，所以我这里就简单地质控一下，由 fastp完成。

质控

#安装fastp
# get source (you can also use browser to download from master or releases)
git clone https://github.com/OpenGene/fastp.git

# build
cd fastp
make

# Install
sudo make install

#没有root权限的话可以使用conda安装或者下载二进制文件
# note: the fastp version in bioconda may be not the latest
conda install -c bioconda fastp

# this binary was compiled on CentOS, and tested on CentOS/Ubuntu
wget http://opengene.org/fastp/fastp
chmod a+x ./fastp

之后批量进行质控，一条龙服务，十分方便了。

#!/bin/bash
for i in {01..45}
do
~/biosoft/fastp/bin/fastp -i /database/Bna-rna-seq-analyis/data/Cole_M838-01-T${i}_good_1.fq.gz -o /database/Bna-rna-seq-analyis/data/T${i}_good_1.fq.gz -I /database/Bna-rna-seq-analyis/data/Cole_M838-01-T${i}_good_2.fq.gz -O /database/Bna-rna-seq-analyis/data/T${i}_good_2.fq.gz
done

比对

比对软件有很多，主流的主要有 HISAT2以及 STAR。

STAR比对

#Download the latest release from and uncompress it
# Get latest STAR source from releases
wget https://github.com/alexdobin/STAR/archive/2.7.2b.tar.gz
tar -xzf 2.7.2b.tar.gz
cd STAR-2.7.2b

# Alternatively, get STAR source using git
git clone https://github.com/alexdobin/STAR.git

#Compile under Linux
# Compile
cd STAR/source
make STAR

STAR使用格式如下：

STAR --option1-name option1-value --option2-name option2-value ...

STAR的文档写的非常详细，这也是一个好软件的标志之一，如果一个软件的文档都写的不详细的话，可以考虑换一个软件。STAR参数非常多，尤其现在还支持三代测序数据。有兴趣的话可以自己去研究 文档。

STAR构建索引

#!/bin/bash

set -e
set -u
set -o pipefail


fasta=/database/reference/Bna_NY7/NY7-chr/NY_V2_Chr.fasta
gtf=/database/reference/Bna_NY7/NY7-chr/NY7_V2_Chr.gtf
STAR=/home/taoyan/biosoft/STAR/bin/Linux_x86_64_static/STAR

# Generate index Of Brassica napus

$STAR --runMode genomeGenerate \
--runThreadN 32 \
--genomeFastaFiles $fasta \
--genomeDir /database/reference/Bna_NY7/NY7-chr/NY7.Star.Index \
--sjdbGTFfile $gtf \
--genomeSAindexNbases 13 \
--sjdbOverhang 149

构建索引常用命令参数如下：

  • –runMode genomeGenerate：运行模式为构建索引
  • –runThreadN：设置线程数
  • –genomeDir：构建的索引文件存放路径，需要提前创建该文件夹
  • –genomeFastaFiles：参考基因组文件
  • –sjdbGTFfile：注释文件gtf，推荐在这一步添加，虽然可以在后续比对的时候添加
  • –genomeSAindexNbases 对于非常小的参考基因组，此参数需要按比例缩小，其计算方式为min(14, log2(GenomeLength)/2 - 1)，所以一般为14，但是小基因组比如1M的话，此值应为9，这里按照我的参考基因组大小，应该为14，但是我运行的时候提示我最优值为13，暂时不清楚为何提示我这个值。此值一般在10~15之间，越大需要越多的内存，相应的话速度也越快。
  • –sjdbOverhang：读段长度-1，默认为99，我的序列是PE150，所以此值应为149，文档显示99可以达到最理想的效果，所以不知道具体的建库方式的话，此参数设为99就行了。

批量比对

#!/bin/bash

set -e
set -u
set -o pipefail


fasta=/database/reference/Bna_NY7/NY7-chr/NY_V2_Chr.fasta
gtf=/database/reference/Bna_NY7/NY7-chr/NY7_V2_Chr.gtf
STAR=/home/taoyan/biosoft/STAR/bin/Linux_x86_64_static/STAR

seq_data=/database/Bna.rnaseq.NY7.analysis/data/seq.data
result=/database/Bna.rnaseq.NY7.analysis/result/Star/Star_map_bam

# sequence mapping use STAR

for i in {01..45}
do
echo " ******* Process Sample T${i} ******* "
$STAR --runThreadN 64 \        #线程数
--genomeDir /database/reference/Bna_NY7/NY7-chr/NY7.Star.Index/ \  #索引目录
--readFilesCommand zcat \   #执行解压缩，输入数据格式为gz文件
--readFilesIn $seq_data/T${i}_good_1.fq.gz \  #输入文件
$seq_data/T${i}_good_2.fq.gz \
--outSAMtype BAM SortedByCoordinate \  #指定输出bam文件并排序
--outFileNamePrefix $result/Bna_Star_T${i}. \  #输出文件路径以及前缀
--outBAMsortingThreadN 20 \   #排序线程数
--quantMode TranscriptomeSAM GeneCounts #将上面基因组比对的信息转化为转录本比对的信息，quantMode TranscriptomeSAM will output alignments translated into trancript coordinates,为了后续使用RSEM定量准备输入文件
done

一些常用参数解释：

  • –runThreadN：线程数
  • –runMode alignReads : 默认就是比对模式，所以这里我没填写
  • –genomeDir：索引文件所在文件夹，就是上面建立索引时的文件夹
  • –readFilesCommand zcat：由于现在大部分序列数据都是压缩形式的，所以需要解压缩，我的数据是gz格式的，所以这里是zcat
  • –readFilesIn：序列数据
  • –outSAMtype BAM SortedByCoordinate：输出文件为bam文件，并且进行了排序，这里应该是按照pos排序的。
  • –outBAMsortingThreadN：SAM文件排序生成BAM文件时的线程数，不要设置太大，不然会报错的。
  • –quantMode TranscriptomeSAM GeneCounts：这个参数需要特别注意，如果你后续要使用RSEM进行定量，就需要添加TranscriptomeSAM，不然生成的BAM文件RSEM无法使用的，GeneCounts生成用于后续featureCounts分析的输入文件

STAR还有其它很多参数，讲不完的，基本来说一般的转录组分析上面这些参数够用了。

由于后续定量我们会用到RSEM，所以这里特别注意需要加参数--quantMode TranscriptomeSAM GeneCounts，这样生成的bam文件可用于后续的RSEM定量以及featureCounts计数。

这里比对的时候我遇到了我问题:

"FATAL ERROR: could not create output file ... Solution: check that you have permission to write this file"

Google搜索之后发现是系统限制了最大文件打开的数目，默认是1024，我改为了10000，后续比对正常运行。

#先查看你的系统打开文件限制数目
$ ulimit -n
1024
#设置一个更大的数目
$ ulimit -n 10000

HISAT2比对

#下载源文件
wget http://ccb.jhu.edu/software/hisat2/downloads/hisat2-2.0.0-beta-source.zip
unzip hisat2-2.0.0-beta-source.zip
cd hisat2-2.0.0-beta-source
make

HISAT2用法如下：

hisat2 [options]* -x <hisat2-idx> {-1 <m1> -2 <m2> | -U <r> | --sra-acc <SRA accession number>} [-S <hit>]

更详细用法见 文档

HISAT2建立索引

我这里遇到一个问题，建立索引之后比对的时候报错：

(ERR): hisat2-align died with signal 11 (SEGV) (core dumped)

Google发现是建立的索引出问题了，HISAT2建立索引的时候会根据基因组大小来建立小索引还是大索引，一般4G以下的基因组都会默认建立小基因组，32位的，以.ht2为后缀的，我这里报错了，所以我特意设置参数建立了64位的以.ht21为后缀的大索引，后续比对就没问题了

#!/bin/bash

set -e
set -u
set -o pipefail


##构建exons以及splice_sites文件，由于我的注释文件是GFF3文件,所以需要转化成GTF

~/biosoft/gffread/gffread NY7_V2_Chr.gff3 -T -o NY7_V2_Chr.gtf
~/biosoft/hisat2-2.0.0-beta/extract_exons.py NY7_V2_Chr.gtf > NY7_V2_Chr.exons
~/biosoft/hisat2-2.0.0-beta/extract_splice_sites.py NY7_V2_Chr.gtf > NY7_V2_Chr.ss

##构建索引,这里需要提供exons，splite_sites文件，参考基因组文件以及输出文件前缀
##Hisat2-build可以构建任意大小的参考基因组索引，小于4G的参考基因组，构建出来的索引以.ht2为后缀的文件，更大的参考基因组是以.ht21为后缀的文件

~/biosoft/hisat2-2.0.0-beta/hisat2-build -p 64 --large-index --ss NY7_V2_Chr.ss --exon NY7_V2_Chr.exon NY_V2_Chr.fasta NY7.Hisat2.Index/NY7_V2_Chr_Hisat2_index 2 > Hisat_build_index.log

批量比对

#!/bin/bash

set -e
set -u
set -o pipefail

hisat2=/home/taoyan/biosoft/hisat2-2.0.0-beta/hisat2
index=/database/reference/Bna_NY7/NY7-chr/NY7.Hisat2.Index
seq_data=/database/Bna.rnaseq.NY7.analysis/data/seq.data
result=/database/Bna.rnaseq.NY7.analysis/result/Hisat2
samtools=/home/taoyan/biosoft/samtools1.9/bin/samtools

for i in {01..45}
do
$hisat2 -t -p 64 --dta -x $index/NY7_V2_Chr_Hisat2_index -1 $seq_data/T${i}_good_1.fq.gz -2 $seq_data/T${i}_good_2.fq.gz | $samtools sort -@ 8 -O bam -o $result/T${i}.sorted.bam
done

具体参数就不详讲了，具体见 文档。

后续主要进行定量了，下次讲。

RNA-seq 比对
Researcher

I am a PhD student of Crop Genetics and Breeding at the Zhejiang University Crop Science Lab. My research interests covers a range of issues:Population Genetics Evolution and Ecotype Divergence Analysis of Oilseed Rape, Genome-wide Association Study (GWAS) of Agronomic Traits.

comments powered by Disqus