1.利用fastqc对模拟测序的序列进行质控分析
1.1 使用art-illumina模拟测序，生成高通量数据（art-illumina模拟测序）
nohup art_illumina -ss HS25 -sam -i GCA_000977205.2_Sc_YJM1338_v1_genomic.fna -p -l 100 -f 45 -m 180 -s 20 -o sequencing1&
nohup art_illumina -ss HS25 -sam -i GCA_000977205.2_Sc_YJM1338_v1_genomic.fna -p -l 100 -f 45 -m 3000 -s 20 -o sequencing2&
1.2 安装fastq及质控分析
下载fastqc http://www.bioinformatics.babraham.ac.uk/projects/download.html#fastqc
解压，进入目录赋予权限
chmod 766 fastqc
echo 'export PATH=$PATH:'`pwd`/ >>~/.bashrc #加入到环境变量
source ~/.bashrc
fastqc主要参数：
-o --outdir FastQC生成的报告文件的储存路径
-t --threads 选择程序运行的线程数
-q --quiet 安静运行模式，一般不选这个选项的时候，程序会实时报告运行的状况
fastqc -t 4 sequencing11.fq sequencing12.fq
fastqc -t 4 sequencing21.fq sequencing22.fq
1.3 总结报告
对sequencing11.fq质量评估后的结果进行分析：
1.从下图中看出模拟测序的质量评估都合格。
sequencing11.fq质量评估总结
2.该表格统计了测序reads长度，总共模拟测序得到2745832条reads，GC含量比为38%.
基本模拟测序数据统计表
3.横轴表示reads位点，纵轴表示质量，蓝线表示此次评估的平均质量。从图中看出平均质量在35左右，最低质量都超过30，所以测序结果很好。
每个位点碱基质量图
4.横轴表示质量，纵轴表示reads条数。大多数reads质量在36左右，所以测序结果很好。
质量分布图
-
横轴表示reads位点，纵轴表示四中碱基的比例。从图中看出“A”‘T’分别占31%左右，“G”“C”分别占19%左右，并且四条线平行，检测合格。
A，T，G，C含量百分比图
-
序列平均GC含量分布图，红线代表实际GC含量，蓝线表示理论GC含量，可以看出在这个fq文件中序列平均GC含量在39%左右。
GC含量分布图
7.sequences duplication是指在测序前建库PCR过程中导致的一些序列扩增次数过多导致的。若重复较高则需要进行处理这些dup。从图中看出重复序列高的次数基本没有。
重复序列统计图
2.测序数据与参考基因组比对
下载并安装bowtie2(https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.5.1/bowtie2-2.3.5.1-sra-linux-x86_64.zip/download)，samtools(https://sourceforge.net/projects/samtools/files/samtools/1.9/samtools-1.9.tar.bz2/download)，下载的samtools是tar.bz2
解压samtools使用下面命令，安装和别的软件一样了。
tar jxvf samtools-1.9.tar.bz2
下载酵母基因组索引文件(ftp://igenome:G3nom3s4u@ussd-ftp.illumina.com/Saccharomyces_cerevisiae/NCBI/build3.1/Saccharomyces_cerevisiae_NCBI_build3.1.tar.gz)，索引文件是为了加快比对。也可以使用bowtie2-build创建索引文件，由于酵母基因组文件不大，所以建索引文件速度还ＯＫ，如果建人类的索引文件还是下载吧，bowtie2建索引太慢了。
#建索引文件，如果已经下载了，则不需要该步
bowtie2-build GCA_000977205.2_Sc_YJM1338_v1_genomic.fna ref/ref
使用bowtie2对上述两套模拟测序的数据进行比对
#下面的ref/ref是索引文件，如果是下载的索引文件则换成相应的文件名
bowtie2 -x ref/ref -U sequencing11.fq,sequencing12.fq -S seq1.sam
bowtie2 -x ref/ref -U sequencing21.fq,sequencing22.fq -S seq2.sam
利用samtools对对比结果进行统计分析
samtools stats seq1.sam --ref-seq GCA_000977205.2_Sc_YJM1338_v1_genomic.fna >gca_stats1.txt
samtools stats seq2.sam --ref-seq GCA_000977205.2_Sc_YJM1338_v1_genomic.fna >gca_stats2.txt
samtools统计结果
序列总数5491664,匹配上的reads有5491662,平均质量36.1
#将sam格式文件转为bam格式，bam是sam的二进制文件，文件小且运算快
samtools view -bS ./seq1.sam > ./seq1.bam
samtools view -bS ./seq2.sam > ./seq2.bam
#对bam格式文件进行排序
samtools sort seq1.bam -o seq1_sort.bam
samtools sort seq2.bam -o seq2_sort.bam
#计算测序深度
samtools depth seq1_sort.bam >seq1_sort_depth.txt
samtools depth seq2_sort.bam >seq2_sort_depth.txt
第一列表示染色体，第二列表示碱基位置，第三列表示每个碱基的测序深度
测序深度
3.序列组装(SOAPdenovo软件)
下载安装SOAPdenovo(https://sourceforge.net/projects/soapdenovo2/files/latest/download)
编写SOAPdenovo的配置文件libc.cfg
max_rd_len=100
[LIB]
avg_ins=180
reverse_seq=0
asm_flags=3
rd_len_cutoff=100
rank=1
pair_num_cutoff=3
map_len=32
q1=../sequencing11.fq
q2=../sequencing12.fq
[LIB]
avg_ins=3000
reverse_seq=1
asm_flags=3
rank=2
pair_num_cutoff=5
map_len=35
q1=../sequencing21.fq
q2=../sequencing22.fq
使用下面命令进行序列组装
nohup SOAPdenovo-63mer all -s lib.cfg -K 30 -o SOAPdenovo_out -p 20 & #服务器上设置20个线程，自己电脑上４个或８个
利用blastn将组装结果与参考基因组进行比对
#创建参考基因组的blast数据库
/opt/BioBuilds/bin/makeblastdb -in GCA_000977205.2_Sc_YJM1338_v1_genomic（genome）.fna -input_type fasta -dbtype nucl -title sc_1338_chr_iv -out sc_1338_chr_iv
#blastn对比
nohup /opt/BioBuilds/bin/blastn -db sc_1338_chr_iv -query ../SOAP/secord/SOAPdenovo_out.contig -out contig_blast_outfmt7 -evalue 1e5 -outfmt 7 -max_target_seqs 1 -num_threads 20 &
nohup /opt/BioBuilds/bin/blastn -db sc_1338_chr_iv -query ../SOAP/secord/SOAPdenovo_out.contig -out contig_blast_outfmt6 -evalue 1e5 -outfmt 6 -max_target_seqs 1 -num_threads 20 &
nohup /opt/BioBuilds/bin/blastn -db sc_1338_chr_iv -query ../SOAP/secord/SOAPdenovo_out.scafSeq -out scafffold_blast_outfmt7 -evalue 1e5 -outfmt 6 -max_target_seqs 1 -num_threads 20 &
从对比结果中提取匹配区间信息，并以参考基因组为基准进行可视化绘图。使用在线工具QUAST，将参考基因组和SOAPdenovo组装生成的 scaffold 文件上传.
网站分析出覆盖度为:91.935,总共176条scaffold
summary
quast图
发现基因组完整，但在scaffold中缺少一段，将该出区域放大
放大后缺少的区域
发现是在cp006396.1号染色体，大致区间在414885－1028442，在基因组中截取该段序列
from Bio import SeqIO
data=SeqIO.parse('GCA_000977205.2_Sc_YJM1338_v1_genomic.fna','fasta')
for seq in data:
if seq.id=='CP006396.1':
with open('seq.fa','w') as f:
f.write(">CP006396.1\n")
f.write(str(seq.seq[414885:1028442]))
将序列上传到ReapeatMasker进行分析
网友评论