基因组信息学实验继续学一点吧

实验二：序列组装

1.art illumina模拟双末端测序

短插入片段和长插入片段库：

./art_illumina -ss HS25 -sam -i ./GCF_000146045.2_R64_genomic.fna -p -l 125 -f 10 -m 200 -s 10 -o ./Sc_paired
./art_illumina -ss HS25 -sam -i ./GCF_000146045.2_R64_genomic.fna -p -l 125 -f 10 -m 2500 -s 50 -o ./Sc_matepair

查看结果：

ll -o --block-size=M ./Sc_paired*
ll -o --block-size=M ./Sc_matepair*

2.bowtie2-build创建索引文件

bowtie2-build GCF_000146045.2_R64_genomic.fna Sc_index

bowtie2-build默认情况下将fasta文件换成index的数据库。
$ bowtie2-build <fasta文件> <要生存的索引文件前缀名>

3.fastqc质控分析

mkdir fastqc_out
fastqc -o ./fastqc_out -f fastq -t 10 /usr/bin/art_bin_MountRainier/Sc_paired1.fq /usr/bin/art_bin_MountRainier/Sc_paired2.fq
fastqc -o ./fastqc_out -f fastq -t 10 /usr/bin/art_bin_MountRainier/Sc_matepair2.fq /usr/bin/art_bin_MountRainier/Sc_matepair1.fq

 

4.bowtie2测序序列与参考序列比对

bowtie2 -x ./Sc_index -1 /usr/bin/art_bin_MountRainier/Sc_paired1.fq,/usr/bin/art_bin_MountRainier/Sc_matepair1.fq -2 /usr/bin/art_bin_MountRainier/Sc_paired2.fq,/usr/bin/art_bin_MountRainier/Sc_matepair2.fq -S ./Sc_2sets.sam -p 10

生成.sam文件

双端数据比对结果： 

 

 第一部分描述的是pair-end模式下的一致比对结果，aligned concordantly就是read1和read2同时合理的比对到了基因组/转录组上。

第二部分，pair-end模式下不一致的比对结果，concordantly 0 times就是read1和read2不能同时合理的同时比对到基因组/转录组上。

第三部分就是对剩余reads（既不能concordantly，也不能discordantly 1 time）的单端模式的比对。

 5.samtools比对结果

用来处理SAM/BAM（SAM的二进制格式，用于压缩空间）格式的比对文件的工具，它能够输入和输出SAM（sequence alignment/map：序列比对）格式的文件，对其进行排序、合并、建立索引等处理。

samtools view -b Sc_2sets.sam >Sc_2sets.bam
#格式转换sam>bam
samtools sort Sc_2sets.bam -o Sc_2sets.sorted.bam
#按照序列名称排序，将结果输出到Sc_2sets.sorted.bam
samtools index Sc_2sets_sorted.bam
#建立索引

 6.做统计分析：

samtools stats ./Sc_2sets_sorted.bam > samtools.stat.stats.out
samtools depth ./Sc_2sets_sorted.bam > samtools.stat.depth.out
samtools flagstat ./Sc_2sets_sorted.bam > samtools.stat.flagstat.out
samtools idxstats ./Sc_2sets_sorted.bam > samtools.stat.idxstats.out

7.解读samtools satas统计结果，利用plot-bamstats对输出结果进行可视化

plot-bamstats -p ./plot-bamstats_out/ ./samtools.stat.stats.out

遇到了error：缺少：gunplot，下载安装后重新运行上述代码

conda install -c bioconda gnuplot -y

8.SOAPdenovo-63mer序列组装及结果分析

nohup SOAPdenovo-63mer all -s lib.cfg -K 31 -o SOAPdenovo_out -p 10 &
#NOHUP  & 放在后台运行

此处因为未安装SOAPdenovo而报错，先下载

git clone https://github.com/aquaskyline/SOAPdenovo2.git
cd SOAPdenovo2
make

因为使用的是Ubuntu 18.04.5而出现了一个报错

自16.10起，gcc就默认开启了pie选项，导致编译出的文件的mime是application/x-sharedlib，一般的文件管理器只认application/x-executable，就没把它当成可执行的文件

修改Makefile中gcc

gcc -fno-pie -no-pie

再make就可以了

 

又因为下列文件的存放路径问题，需修改配置文档

 文件路径为

/usr/bin/art_bin_MountRainier

改了一堆问题想必是可以运行了试试吧，就不在后台运行了

SOAPdenovo2/SOAPdenovo-63mer all -s lib.cfg -K 31 -o SOAPdenovo_out -p 1

 不错果然跑完了，小辣鸡真不容易

9.quast将组装结果中的contigs和scaffolds序列的文档分别与参考基因组进行对比

quast-5.0.2/quast.py -o quast_out -r GCF_000146045.2_R64_genomic.fna -g GCF_000146045.2_R64_genomic.gff SOAPdenovo_out.contig
#contigs评估
quast-5.0.2/quast.py -o quast_out -r GCF_000146045.2_R64_genomic.fna -g GCF_000146045.2_R64_genomic.gff SOAPdenovo_out.scafSeq
#scaffolds评估

 10.结果查看就留到有空吧

呜呜呜下班了