@Billy-The-Crescent 2019-06-20T23:11:13.000000Z 字数 3358 阅读 460

计生项目实践-hqj

计算生物学 项目 Protocol 胡泉军

基因组组装

GenomeAssembly：将测序得到的序列片段组装成一个全长的基因组序列

reads：就是我们测序产生的短读序列，通常一代和三代的reads读长在几千到几万bp之间，二代的相对较短，平均是几十到几百bp。

contig：中文叫做重叠群，就是不同reads之间的overlap交叠区，拼接成的序列就是contig

scaffold:这是比contig还要长的序列，获得contig之后还需要构建paired-end或者mate-pair库，从而获得一定片段的两端序列，这些序列可以确定contig的顺序关系和位置关系，最后contig按照一定顺序和方向组成scaffold，其中形成scaffold过程中还需要填补contig之间的空缺。

分为三步：

将reads组装成contig

将contig组装成scaffold

将scaffold组装成基因组

reads to contig

使用方法：基于De Brujin Graph(DBG)的方法

使用软件：Velvet，manual
输入：正向和反向测序的reads序列文件 (fastq格式)
input: SRR292770_1.fastq.gz SRR292770_2.fastq.gz分别代表正向和反向测序数据。

加入velvet的环境变量

export PATH=/data/apps/velvet/velvet_1.2.10:$PATH

首先，指定k-mer的长度，生成k-mer，这里默认使用35。因此，组装命令如下：

velveth out_data_35 35 -fastq.gz -shortPaired -separate SRR292770_1.fastq.gz SRR292770_2.fastq.gz

其中，-shortPaired表示使用的是短的双末端序列，-seperate表示正向序列和反向测序序列是分开存放的。

生成的结果存放在新建的out_data_35目录下。
3. 使用k-mer构建都布灵图：

velvetg out_data_35 -clean yes -exp_cov 21 -cov_cutoff 2.81 -min_contig_lgth 200

最小contig长度设置为200

k-‐mer = 35, expected coverage = 20, coverage cutoff of 2.81

从结果中提取contig的.fasta文件，并删除中间文件

cp out_data_35/contigs.fa SRR292770_unordered.fasta
rm -r out_data_35

优化velvet的de novo组装

使用软件： VelvetOptimiser

加入环境变量，PERL5LIB的路径VelvetOptimiser软件的路径，并确保velvet的路径在环境变量中。

export PATH=/data/apps/velvet/velvet_1.2.10:$PATH
export PERL5LIB=/data/apps/cpan/lib/perl5:/data/apps/cpan/lib64/perl5:/data/apps/cpan/share/perl5
export PATH=/data/apps/velvet/velvet_1.2.10/contrib/VelvetOptimiser-2.2.4:$PATH

运行VelvetOptimiser

VelvetOptimiser.pl -s 33 -e 41 -f '-fastq.gz -shortPaired SRR292770_1.fastq.gz SRR292770_2.fastq.gz' -o '-min_contig_lgth 200' -p SRR292770

VelvelOptimiser尝试了长度在33到41之间的所有奇数的k-mer，并进行优化，得到的结果在SRR292770_data_35
，SRR292770是前缀，和之前的结果进行区分。

将得到的unordered contig进行排序

通过reference genome的序列，将无序的contig序列比对在参考基因组上，从而进行排序

使用软件：abacas

加入环境变量：

export PATH=/data/apps/MUMmer3.23:$PATH
export PERL5LIB=/data/apps/MUMmer3.23/scripts:$PERL5LIB
export PATH=/data/apps/abacas/1.3.1/:$PATH

将所有需要用到的软件的路径都添加到环境变量中。
2. 获取reference genome：

cp /var/www/html/p/2019-ComBioPra/Genome/data/NC_011748.fasta NC_011748.fasta

将参考基因组复制到目录下
3. 进行对比，排序：

abacas.pl -r NC_011748.fasta -q SRR292770_unordered.fasta -p 'nucmer' -c -m -b -o SRR292770

使用abacas.pl脚本进行排序，得到的fasta文件的前缀为SRR292770。

使用Artemis Comparison Tool (ACT)进行contig可视化

软件：ACT
首先将拼接好的contig和reference genome进行blast比对

export PATH=/data/apps/blast/ncbi-blast-2.7.1p/bin:$PATH
makeblastdb -in NC_011748.fasta -dbtype nucl
blastn -query SRR292770.fasta -db NC_011748.fasta -evalue 1 -task megablast -outfmt 6 > SRR292770_NC_011748.crunch

然后将比对好的.crunch文件和fasta文件输入到ACT软件中就可以得出结果。

转录组分析

设置转录组分析需要用到的所有环境变量：

time=/usr/bin/time
export TRINITY_HOME=/data/apps/trinity/trinityrnaseq-Trinity-v2.8.4
bowtieDir=/data/apps/bowtie/bowtie2-2.3.2/bin
samtoolsDir=/data/apps/samtools/1.9/bin
tophatDir=/data/apps/tophat/tophat-2.1.1.Linux_x86_64
gmapDir=/data/apps/gmap/20190315/bin
cufflinksDir=/data/apps/cufflinks/2.2.1
RDir=/data/apps/R/3.5.1/bin
rsemDir=/data/apps/rsem/1.3.1/bin
modulesDir=/data/apps/modules/3.2.10/Modules/3.2.10/bin
javaDir=/data/apps/jdk/jdk1.8.0_131/bin
jellyfishDir=/data/apps/jellyfish/2.2.9
salmonDir=/data/apps/salmon/0.9.1/bin
bowtie1Dir=/data/apps/bowtie/1.2.2
resmDir2=/data/apps/rsem/1.3.1/usr/local/bin
export PATH=$TRINITY_HOME:$bowtieDir:$samtoolsDir:$tophatDir:$gmapDir:$cufflinksDir:$RDir:$rsemDir:$modulesDir:$javaDir:$jellyfishDir:$salmonDir:$PATH
export PATH=$bowtie1Dir:$resmDir2:$PATH
export PATH=/data/apps/R/3.5.1/bin:$PATH
export PATH=/data/apps/blast/ncbi-blast-2.7.1p/bin:$PATH