多态性核SSR的鉴定

---

多态性核SSR的鉴定

所属目录：紫菜作者：星云<XingYun>创建时间：2024/10/28更新时间：2024/11/6URL：

前言

本文主要记录了本人如何使用bwa对测序数据进行mapping，再使用SOAPdenovo2对核序列进行从头组装成scaffolds，再使用CandiSSR寻找多态性核SSR。最后整理成一个文件。

参考文献：Intraspecific phylogeny and genomic resources development for an important medical plant Dioscorea nipponica, based on low-coverage whole genome sequencing data

根据文献中的所用方法步骤来完成：

整理形成的文件表格：

最后所整理的文件（下载Table S4文件查看类似的表格）：

一、使用bwa对测序数据进行mapping

主要是使用bwa软件将测序数据mapping到参考基因组上，以排除线粒体和叶绿体reads,得到核序列。
参考文章：
【比对软件】BWA使用
BWA使用详解

使用示例：

# 假设参考基因组文件为yezoensis.genome.fa，cleaned reads文件为ST1-3_clean_R1.fq.gz和ST1-3_clean_R2.fq.gzref=yezoensis.genome.fa
reads1=ST1-3_clean_R1.fq.gz
reads2=ST1-3_clean_R2.fq.gz# 建立BWA索引
bwa index -a bwtsw $ref# 对clean_R1.fq.gz和clean_R2.fq.gz进行单端序列的mates拼接
bwa mem -t 4 $ref $reads1 $reads2 > ST1-3_aligned.sam # 将SAM格式的文件转换为BAM格式的文件
samtools sort -o ST1-3_aligned.bam ST1-3_aligned.sam# 使用samtools对BAM文件进行索引
samtools index ST1-3_aligned.bam

二、使用SOAPdenovo2对核序列进行从头组装成scaffolds

进行完上一步mapping后，得到了bam文件。之后就是使用SOAPdenovo2软件对核序列进行从头组装成scaffolds。
参考文章：基因组组装—SOAPdenovo2的使用

SOAPdenovo2软件包：
链接：https://pan.baidu.com/s/1IK02W2UEa8v9wGI0BEhsWw
提取码：8uvl

首先进入SOAPdenovo2软件包的目录，再使用命令

SOAPdenovo2的使用需要自己构建配置文件，可以根据软件目录中的示例配置文件"example.config"进行构建
以上面得到的一个bam文件：ZD-1_aligned.bam为例，在当前目录下创建一个config_file文件：ZD-1.config

vim ZD-1.config

添加以下内容（仅供示例，请结合软件使用手册根据自己序列的具体情况进行修改）

max_rd_len=150
[LIB]
avg_ins=350
reverse_seq=0
asm_flags=3
rd_len_cutoff=150
rank=1
pair_num_cutoff=3
map_len=32
b=../ZD-1_aligned.bam

然后直接一站式运行

./SOAPdenovo-63mer all -s ZD-1.config -K 63 -R -o ZD-1

运行的结果文件中，其中有下面两个主要的组装结果文件

*.contig  # contig序列文件
*.scafSeq # scaffold序列文件

我这里要的是scaffold序列文件，可以把它们移动到新的文件夹中，方便查看和使用

mkdir scafSeq
mv *.scafSeq scafSeq

至于contig序列文件，我们也可以把它们移动到新的文件夹中。

mkdir contig
mv *.contig contig

三、使用CandiSSR寻找多态性核SSR

CandiSSR的项目地址：CandiSSR

3.1. 安装CandiSSR软件的准备

根据项目的readme文档，在安装此软件前，需要先保证环境的配置，也就是安装所要的依赖软件(也可以先安装此CandiSSR项目，再使用项目中提供的软件去进行配置)

安装完所要的依赖软件后，就可以启动setup.sh，将所依赖的软件路径配置好后，形成CandiSSR.pl脚本。

CandiSSR.pl脚本文件内容（我们需要做的就是将软件的可执行路径添加进去，一般perl软件在实验室服务器中都有安装，也就不必要额外再安装了）：

补充，在我安装Primer3软件时，使用项目自带的软件编译时报错，无法形成执行文件primer3_core。所以我就从github上下载了primer3软件，能正常使用。

GIThub上的Primer3下载地址

注意：如果在后面运行时有报错或者没有形成最后的文件。可以检查一下CandiSSR.pl脚本中是不是正确配置好了所需软件路径。总之就是结合报错信息来对脚本内容进行修改。

3.2. 运行CandiSSR时的准备

形成正确的CandiSSR.pl文件后，就可以使用命令来运行了，按照软件的项目文档介绍来运行

先形成一个配置文件 .ctl,(第一行为参考序列，其余行为要分析鉴定的序列)

待分析鉴定的序列为上一步形成的scaffold序列文件，又因为按照示例的配置文件中，序列文件都都是以fasta为后缀，所以我们可以将scaffold序列文件转换为fasta格式(改一下后缀名就行，也许可以不改，因为原来的格式就是fasta格式，只是后缀不一样)

运行示例（使用默认参数）

perl CandiSSR.pl -i ref_haitanensis.ctl -o haitanensis -p haitanensis 2>&1 | tee haitanensis_output.log

运行完后，一般会得到三个文件

3.3. 整理得到的结果文件

也就是将_PolySSRs.txt和_Designed_Primers.txt文件合并在一个表格中

四、统计Contig的数量和N50值

前面第二步：使用SOAPdenovo2对核序列进行从头组装成scaffolds。得到了contig序列文件。而现在这一步则是使用前面得到的contig序列文件，统计得到contig的数量和N50值。

参考文章：20220518_基因组contig与scaffold的N50大小统计

直接使用python程序来计算：

def calculate_n50(contig_lengths):total_length = sum(contig_lengths)contig_lengths.sort(reverse=True)n50 = 0for length in contig_lengths:n50 += lengthif n50 >= total_length / 2:return lengthdef parse_fasta(file_path):contig_lengths = []with open(file_path, 'r') as file:current_length = 0for line in file:if line.startswith('>'):if current_length > 0:contig_lengths.append(current_length)current_length = 0else:current_length += len(line.strip())if current_length > 0:contig_lengths.append(current_length)return contig_lengths# 使用方法
file_path = r''  # 替换为你的FASTA文件路径
contig_lengths = parse_fasta(file_path)
n50 = calculate_n50(contig_lengths)
num_contigs = len(contig_lengths)print(f'Contig数量: {num_contigs}')
print(f'N50长度: {n50}')

或者可以直接使用一个python库：assembly_stats来统计contig文件的相关信息（如果要统计scaffold的相关信息则要换scaffold文件）

#先安装python库
pip install assembly_stats#在终端运行
assembly_stats 序列文件路径

五、使用OmicStudio进行数据可视化

暂无（因为主要是对多态性核SSR进行鉴定，得到文件就行，对于可视化也没有要求，有空再看看）

六、总结

本文主要学会的是如何对测序数据进行多态性核SSR的鉴定。
从第一步的使用bwa软件将测序数据映射到参考序列中，获得仅包含核reads的bam数据；
到第二步的使用SOAPdenovo2程序将bam数据进行从头组装成scaffolds；
再到最后的使用CandiSSR软件使用默认参数来识别候选多态性SSRs，这一步涉及了许多依赖软件的安装，比较麻烦。

另外，也学会统计了contig的数量和N50值。

2024/11/6