nature新冠病毒基因序列测定流程

写在前头

新冠肺炎导致了严重的全球健康危机,目前针对新冠病毒的检测,预防手段都已经成熟,以生物学角度来看,新冠病毒通过将其遗传物质(RNA)注射进入人体细胞,从而完成自身的繁殖。

本教程将尝试利用sixoclock平台在新窗口打开复现2020年武汉发表的揭示新冠病毒遗传物质(RNA序列)的nature文章《A new coronavirus associated with human respiratory disease in China》在新窗口打开(以下简称新冠nature文章)中的生物信息学流程,

并比较以常规方法在linux系统上搭建流程的分析策略与使用sixoclock平台的分析策略有什么异同,进一步体会云计算协作技术如何加速生物医疗数据分析的速度,降低使用门槛,提高研究人员的工作效率,使分析更专注于问题本身而不是工具等的安装与使用。

通过本案例,你将了解:

  • 一个典型的生物信息过程是什么?
  • sixoclock平台提供怎样的服务
  • sixoclock如何提高数据分析效率

新冠nature文章简介

本次节选的新冠nature文章在新窗口打开是发表在nature的一篇揭示新冠病毒遗传物质的研究性文章,文章通过对新冠病人进行基因测序,利用生物信息分析方法,确定了新冠病毒的核酸序列,从而为后续针对新冠病毒的筛查,预防与治疗奠定了基础。

文章分析流程

img

可以看到,文章使用了较多工具,为了复现文章,我们先列出所有用到的算法/软件(见清单一),与待分析数据(见清单二)。

清单一——软件清单

软件版本网站
Trimmomatic0.39http://www.usadellab.org/cms/?page=trimmomatic
Bowtie22.2.9http://bowtie-bio.sourceforge.net/bowtie2/index.shtml
Megahitv.1.1.3https://github.com/voutcn/megahit
Trinityv2.8.5https://github.com/trinityrnaseq/trinityrnaseq
rsem1.3.3https://github.com/deweylab/RSEM/
BLASTn2.10.0https://blast.ncbi.nlm.nih.gov/Blast.cgi
Diamond2.0.6https://github.com/bbuchfink/diamond
Samtools1.11http://www.htslib.org/doc/samtools.html
PhyMLv3.0http://github.com/stephaneguindon/phyml

清单二 ——流程输入数据

数据类型下载链接
hg38.fa.gz人参考基因组http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz
viral.1.1.genomic.fna.gz已知各种病毒的基因序列https://ftp.ncbi.nlm.nih.gov/refseq/release/viral/viral.1.1.genomic.fna.gz
viral.1.protein.faa.gz已知各种病毒的蛋白序列http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz
WH_R1.fastq.gz新冠病人基因测序数据https://ftp.cngb.org/pub/Dataset/datadis19/SRP245409/SRR10971381/WH_R1.fastq.gz
WH_R2.fastq.gz新冠病人基因测序数据https://ftp.cngb.org/pub/Dataset/datadis19/SRP245409/SRR10971381/WH_R2.fastq.gz
TruSeq3-PE.fa测序设备标识序列(测序接头)http://www.sixoclock.net/resources/data/NGS/Homo_sapiens/Others/adapters/TruSeq3-PE.fa

基于sixoclock平台复现文章

sixoclock平台的软件仓库目前托管了大量数据分析算法/软件,依托于已有工作,我们可以轻松获取清单一中的软件,

登陆sixoclock平台,查找软件Trimmomatic, 并通过可视化界面完成软件参数设置(详情可参考六点了入门教程)

image-20210620184153507

在线完成该软件参数配置,

image-20210620184312036

其它软件类似,sixoclock平台的软件基于docker,因而可以直接在安装有docker/singlarity/udocker 的机器上运行,可省去用户繁琐的软件安装与编译过程,实现一键下载,一键运行。

这里,我们整理了一个完整的文章复现流程,最终复现的流程结构如图:

image-20210620184312036

您可前往新冠nature复现流程在新窗口打开查看数据分析流程结构,点击设置运行按钮以尝试通过可视化点击来替换参数,进而实现个性化分析,最后下载流程与配置文件,本地运行以利用已有工作快速完成文章的复现。

当然,您也可以参考该流程,实现您自有数据的分析,充分享受sixoclock协作平台带给你的便利性。

详细过程如下:

首先,

下载流程描述文件在新窗口打开coronavirus_re-emerged_pipeline.cwl和参数配置文件coronavirus_re-emerged_pipeline.job.yaml,可下载我们准备好的测试demo在新窗口打开

或者按如下设置流程参数文件(保存为coronavirus_re-emerged_pipeline.yml):

threads: 2  # type "int" (optional)
out_pre_1: blastn_megahitfa.out  # type "string" (optional)
outfmt_3: 6  # type "int" (optional)
evalue_3: 0.0000000001  # type "float" (optional)
outfmt_2: 6  # type "int" (optional)
out_pre: blastn_trinityfa.out  # type "string" (optional)
evalue_2: 0.0000000001  # type "float" (optional)
o_1: blastx_megahitfa.out  # type "string" (optional)
o: blastx_trinityfa.out  # type "string" (optional)
evalue_1: 0.00001  # type "float" (optional)
outfmt_1: 6  # type "int" (optional)
outfmt: 6  # type "int" (optional)
evalue: 0.00001  # type "float" (optional)
max_memory: 10G  # type "string"
bt2_index_base: hostfa  # type "string" (optional)
un_conc: rmhost.clean.fq.gz  # type "string" (optional)
end_mode: PE  # type "string"
illuminaclip: '2:30:10'  # type "string"
infile:  # type "File"
    class: File
    path: http://www.sixoclock.net/resources/data/NGS/SARS-COV-2/Panel/viral.1.protein.faa
input_file:  # type "File"
    class: File
    path: http://www.sixoclock.net/resources/data/NGS/SARS-COV-2/Panel/viral.1.1.genomic.fna
reference_in:  # type "File"
    class: File
    path: http://www.sixoclock.net/resources/data/NGS/Homo_sapiens/Reference/hg19.chrM.fasta
input_read2_fastq_file:  # type "File" (optional)
    class: File
    path: http://www.sixoclock.net/resources/data/NGS/SARS-COV-2/RNA_Seq/nCoVR.WH_2.fastq.gz
input_adapters_file:  # type "File"
    class: File
    path: http://www.sixoclock.net/resources/data/NGS/Homo_sapiens/Others/adapters/TruSeq3-PE.fa
input_read1_fastq_file:  # type "File"
    class: File
    path: http://www.sixoclock.net/resources/data/NGS/SARS-COV-2/RNA_Seq/nCoVR.WH_1.fastq.gz
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

然后,

在linux操作系统,安装sixbox ,并执行如下命令运行流程:

sixbox run coronavirus_re-emerged_pipeline.cwl coronavirus_re-emerged_pipeline.job.yaml
1

上述命令将会依次拉取所需的docker镜像,并按照配置文件开始执行命令,进行数据分析过程。

常规方法复现流程

常规分析方法进行生物信息分析过程,一般包含软件安装编译demo测试命令脚本拟写参数设置流程组合几部分。

首先,我们依次下载安装清单一所列软件,鉴于软件过多,此处仅列出几个代表性的软件

软件安装与编译

安装Trimmomatic

TRIMMOMATIC_VERSION=0.39
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-${TRIMMOMATIC_VERSION}.zip && \
    unzip Trimmomatic-${TRIMMOMATIC_VERSION}.zip
1
2
3

安装Trinity

 apt-get -qq update && apt-get -qq -y install \
    automake \
    build-essential \
    bzip2 \
    cmake \
    curl \
    default-jre \
    fort77 \
    ftp \
    g++ \
    gcc \
    gfortran \
    git \
    libblas-dev \
    libbz2-dev \
    libcairo2-dev \
    libcurl4-openssl-dev \
    libdb-dev \
    libghc-zlib-dev \
    libjpeg-dev \
    liblzma-dev \
    libncurses-dev \
    libncurses5-dev \
    libpcre3-dev \
    libpng-dev \
    libreadline-dev \
    libreadline-dev \
    libssl-dev \
    libtbb-dev \
    libx11-dev \
    libxml2-dev \
    libxt-dev \
    libzmq3-dev \
    make \
    nano \
    perl \
    pkg-config \
    python3 \
    python3-dev \
    python3-distutils \
    python3-pip \
    python3-setuptools \
    rsync \
    texlive-latex-base \
    tzdata \
    unzip \
    wget \
    x11-common \
    zlib1g-dev

## Perl stuff
curl -L https://cpanmin.us | perl - App::cpanminus
cpanm install DB_File
cpanm install URI::Escape

## set up tool config and deployment area:
SRC=/usr/local/src
BIN=/usr/local/bin

#####
# Install R
R_VERSION=R-3.6.3

curl https://cran.r-project.org/src/base/R-3/$R_VERSION.tar.gz -o $R_VERSION.tar.gz && \
    tar xvf $R_VERSION.tar.gz && \
    cd $R_VERSION && \
	./configure && make && make install

R -e 'install.packages("BiocManager", repos="http://cran.us.r-project.org")'
R -e 'BiocManager::install("tidyverse")'
R -e 'BiocManager::install("edgeR")'
R -e 'BiocManager::install("DESeq2")'
R -e 'BiocManager::install("ape")'
R -e 'BiocManager::install("ctc")'
R -e 'BiocManager::install("gplots")'
R -e 'BiocManager::install("Biobase")'
R -e 'BiocManager::install("qvalue")'
R -e 'BiocManager::install("goseq")'
R -e 'BiocManager::install("Glimma")'
R -e 'BiocManager::install("ROTS")'
R -e 'BiocManager::install("GOplot")'
R -e 'BiocManager::install("argparse")'
R -e 'BiocManager::install("fastcluster")'
R -e 'BiocManager::install("DEXSeq")'
R -e 'BiocManager::install("tximport")'
R -e 'BiocManager::install("tximportData")'

LD_LIBRARY_PATH=/usr/local/lib

apt-get install -y cython3
## Python 3 stuff
ln -sf /usr/bin/python3 /usr/bin/python

## some python modules
pip3 install numpy
pip3 install git+https://github.com/ewels/MultiQC.git
pip3 install HTSeq

## bowtie
cd $SRC
wget https://sourceforge.net/projects/bowtie-bio/files/bowtie/1.2.1.1/bowtie-1.2.1.1-linux-x86_64.zip/download -O bowtie-1.2.1.1-linux-x86_64.zip && \
    unzip bowtie-1.2.1.1-linux-x86_64.zip && \
	mv bowtie-1.2.1.1/bowtie* $BIN

## RSEM
mkdir /usr/local/lib/site_perl
cd $SRC
wget https://github.com/deweylab/RSEM/archive/v1.3.3.tar.gz && \
    tar xvf v1.3.3.tar.gz && \
    cd RSEM-1.3.3 && \
    make && \
    cp rsem-* convert-sam-for-rsem $BIN && \
    cp rsem_perl_utils.pm /usr/local/lib/site_perl/ && \
    cd ../ && rm -r RSEM-1.3.3

## Kallisto
cd $SRC
wget https://github.com/pachterlab/kallisto/releases/download/v0.46.1/kallisto_linux-v0.46.1.tar.gz && \
    tar xvf kallisto_linux-v0.46.1.tar.gz && \
    mv kallisto/kallisto $BIN

## FASTQC
cd $SRC
wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.5.zip && \
    unzip fastqc_v0.11.5.zip && \
    chmod 755 /usr/local/src/FastQC/fastqc && \
    ln -s /usr/local/src/FastQC/fastqc $BIN/.

# blast
cd $SRC
wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.5.0/ncbi-blast-2.5.0+-x64-linux.tar.gz && \
    tar xvf ncbi-blast-2.5.0+-x64-linux.tar.gz && \
    cp ncbi-blast-2.5.0+/bin/* $BIN && \
    rm -r ncbi-blast-2.5.0+

## Bowtie2
cd $SRC
wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.4.1/bowtie2-2.3.4.1-linux-x86_64.zip/download -O bowtie2-2.3.4.1-linux-x86_64.zip && \
    unzip bowtie2-2.3.4.1-linux-x86_64.zip && \
    mv bowtie2-2.3.4.1-linux-x86_64/bowtie2* $BIN && \
    rm *.zip && \
    rm -r bowtie2-2.3.4.1-linux-x86_64

## Samtools
wget https://github.com/samtools/samtools/releases/download/1.10/samtools-1.10.tar.bz2 && \
    tar xvf samtools-1.10.tar.bz2 && \
    cd samtools-1.10 && \
    ./configure && make && make install

## Jellyfish
wget https://github.com/gmarcais/Jellyfish/releases/download/v2.2.7/jellyfish-2.2.7.tar.gz && \
    tar xvf jellyfish-2.2.7.tar.gz && \
    cd jellyfish-2.2.7/ && \
    ./configure && make && make install

## Salmon
cd $SRC
SALMON_VERSION=1.4.0
wget https://github.com/COMBINE-lab/salmon/releases/download/v${SALMON_VERSION}/Salmon-${SALMON_VERSION}_linux_x86_64.tar.gz && \
    tar xvf Salmon-${SALMON_VERSION}_linux_x86_64.tar.gz && \
    ln -s $SRC/Salmon-latest_linux_x86_64/bin/salmon $BIN/.

## FeatureCounts
wget https://sourceforge.net/projects/subread/files/subread-2.0.0/subread-2.0.0-Linux-x86_64.tar.gz/download -O subread-2.0.0-Linux-x86_64.tar.gz && \
    tar xvf subread-2.0.0-Linux-x86_64.tar.gz && \
    cp -r subread-2.0.0-Linux-x86_64/bin/* $BIN/

## Hisat2
wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/downloads/hisat2-2.1.0-Linux_x86_64.zip && \
    unzip hisat2-2.1.0-Linux_x86_64.zip && \
    cp hisat2-2.1.0/hisat2* $BIN/

## GMAP
GMAP_VERSION=2017-11-15
cd $SRC
GMAP_URL="http://research-pub.gene.com/gmap/src/gmap-gsnap-$GMAP_VERSION.tar.gz" && \
    wget $GMAP_URL && \
    tar xvf gmap-gsnap-$GMAP_VERSION.tar.gz && \
    cd gmap-$GMAP_VERSION && ./configure --prefix=`pwd` && make && make install && \
    cp bin/* $BIN/

# patch for salmon
cd $SRC
ln -sf /usr/local/src/salmon-latest_linux_x86_64/bin/salmon $BIN/.

# blat
RUN wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/blat/blat -P $BIN && \
    chmod a+x $BIN/blat

## Picard tools
cd $SRC
wget https://github.com/broadinstitute/picard/releases/download/2.23.3/picard.jar
PICARD_HOME=$SRC

####
## GATK4 installation
cd $SRC
GATK_VERSION=4.1.8.0
wget https://github.com/broadinstitute/gatk/releases/download/${GATK_VERSION}/gatk-${GATK_VERSION}.zip && \
    unzip gatk-${GATK_VERSION}.zip

GATK_HOME=$SRC/gatk-${GATK_VERSION}

## STAR
STAR_VERSION=2.7.8a
STAR_URL="https://github.com/alexdobin/STAR/archive/${STAR_VERSION}.tar.gz" &&\
    wget -P $SRC $STAR_URL &&\
        tar -xvf $SRC/${STAR_VERSION}.tar.gz -C $SRC && \
            mv $SRC/STAR-${STAR_VERSION}/bin/Linux_x86_64_static/STAR /usr/local/bin

## Trinity
cd $SRC
TRINITY_VERSION="2.12.0"
TRINITY_CO="a1431ffa1d87236dba814ddf0bb2e5c1f45c43d2"

cd $SRC

git clone --recursive https://github.com/trinityrnaseq/trinityrnaseq.git && \
    cd trinityrnaseq && \
    git checkout ${TRINITY_CO} && \
    git submodule init && git submodule update && \
    git submodule foreach --recursive git submodule init && \
    git submodule foreach --recursive git submodule update && \
    rm -rf ./trinity_ext_sample_data && \
    make && make plugins && \
    make install && \
    cd ../ && rm -r trinityrnaseq

TRINITY_HOME=/usr/local/bin

PATH=${TRINITY_HOME}:${PATH}

# some cleanup
cd $SRC
rm -r ${R_VERSION} *.tar.gz *.zip *.bz2

apt-get -qq -y remove git && \
    apt-get -qq -y autoremove && \
    apt-get clean && \
    rm -rf /var/lib/apt/lists/* /var/log/dpkg.log

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241

其它软件安的装方法大致相同,可以看到,部分软件依赖复杂,仅安装软件这一步就耗费大量时间,且,对于一般用户,摸索在自己系统环境中正确安装软件是一件较为困难的事情。

流程搭建

makeblastdb \
 -dbtype nucl \
 -in WTX.gene.fa \
 -parse_seqids \
 -title WTX.gene \
 -input_type fasta \
 -out WTX.gene


bowtie2-build \
 --threads 2 \
 hg19.chrM.fasta
 hg19.chrM

diamond \
 makedb \
 --in WTX.database.fa
 --db WTX.database
 
java -jar /home/6oclock/bin/trimmomatic-0.39.jar
 PE \
 -threads 2 \
 -phred64 \
 NA12878.WGS_1.fastq.gz \
 NA12878.WGS_2.fastq.gz \
 NA12878.WGS_1.fastq.trimmed.fastq \
 NA12878.WGS_1.fastq.trimmed.unpaired.fastq \
 NA12878.WGS_2.fastq.trimmed.fastq \
 NA12878.WGS_2.fastq.trimmed.unpaired.fastq \
 ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
 HEADCROP:34 \
 LEADING:3 \
 TRAILING:3
 
bowtie2
 --dpad 0 \
 --gbar 99999999 \
 --no-discordant \
 --no-mixed \
 --sensitive \
 -p 4 \
 --un-conc-gz unmap.fq.gz \
 -x hg19.chrX.fasta_bowtie2_index/hg19.chrM \
 -S .sam \
 -1 \
 NA12878.WGS-100K_1.fastq.gz \
 -2 \
 > .bowtie2_stdout
 
megahit \
 --out-prefix final \
 -1 nCoVR.WH-100K_1.fastq.gz \
 -2 nCoVR.WH-100K_2.fastq.gz \
 --k-list 21,29,39,59,79,99,119,141 \
 --min-count 2 \
 -t 2 \
 -o ./megahit_out
 
 Trinity \
 --CPU 4 \
 --seqType fq \
 --left NA12878.WGS_1.fastq \
 --max_memory 10G \
 --right NA12878.WGS_2.fastq \
 --output trinity
 
rsem-prepare-reference \
 --gtf hg19.chrM.refGene.gtf \
 -p 1 \
 hg19.chrM.fasta \
 hg19.chrM \
 > rsem-prepare-reference.log

rsem-calculate-expression \
 --alignments \
 --paired-end \
 nCoVR.WH_bowtie2.fastq.sam \
 nCoVR.contig.rsem_index/rsemRNA \
 test \
 > test.rsem.stdout.log

diamond blastx \
 -d WTX.database.dmnd \
 --threads 1 \
 --outfmt 6 \
 -q WTX.gene.fa \
 --out ncov.diamondblastx.out \
 --evalue 0.00001
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88

依次替换上述代码中的输入数据,运行shell脚本完成流程运行。