重测序数据处理得到vcf文件

---

重测序数据处理

所属目录：紫菜创建时间：2024/7/8作者：星云<XingYun>更新时间：2024/7/22URL：https://blog.csdn.net/2301_78630677/article/details/140570255

前言

本文主要是记录了重测序数据处理流程：对于送往公司的DNA样品进行重测序后，得到的数据进行一系列处理得到vcf文件。

有关于重测序数据处理的理论以及其它相关了解内容，若想了解，请查看另一篇为此专门写的博客：补充——重测序数据处理的理论以及其它相关了解内容

数据处理程序流程图

！！！注意： 后续处理流程的代码仅供参考，请根据自己的实际情况自行修改！

1. 数据是rawdata，需用fastp对数据进行质控和过滤

有关使用fastp对数据进行质控和过滤，推荐阅读：fastp：数据质控 + 过滤

fastp -i GQS10_1_T_1.fq.gz -I GQS10_1_T_2.fq.gz -o GQS10_1_T_1.clean.fq.gz -O GQS10_1_T_2.clean.fq.gz

2. 利用getorganelle软件组装叶绿体基因组

有关使用getorganelle软件组装叶绿体基因组，推荐阅读：零基础教程 | 叶绿体基因组组装 - GetOrganelle

get_organelle_from_reads.py -1 GQS10_1_T_1.clean.fq.gz -2 GQS10_1_T_2.clean.fq.gz
-F embplant_pt -o GQS10_1_T-cp -R 10 -t 30 -k 21,45,65,85,105

3. 检查基因组大小，确认是否完整，然后和已知的红毛菜科叶绿体基因组一起构树

有关如何构树，请参考之前我的另一篇博客：实验篇——多序列比对，构树 , 只不过用的是muscle比对，但具体流程差不多。

有关于mafft比对与muscle比对。推荐查看这篇文章：4-【mafft、muscle】的安装和使用

先用mafft比对
trimal修剪改格式
iqtree2构树

4. 根据树形结果挑选坛紫菜个体，为了后续分析方便，对所有个体进行重命名

在这一步，标记好哪些文件是坛紫菜，用于后续的挑选（如果此时就将需要的样品直接挑选出来进行后续的分析，那么第11步就直接使用第一个脚本文件；如果此时只是记录，没有直接挑选出来，而是一起进行后续分析，那么第11步使用第二个脚本文件）

5. 给参考基因组建索引

T2T_genome.fa为参考基因组；
这段代码是用于对一个名为 T2T_genome.fa 的基因组序列文件进行索引和字典创建。

bwa index T2T_genome.fa  
samtools faidx T2T_genome.fa
samtools dict T2T_genome.fa>T2T_genome.dict

6. clean reads 比对到参考基因组生成bam文件并对bam文件进行排序

这段代码的主要目的是使用bwa对clean.fq.gz文件（_1.clean.fq.gz与_2.clean.fq.gz）进行比对，并将结果保存为bam格式，然后进行排序

工作路径：在存放.clean.fq.gz数据的目录

for fname in *_1.clean.fq.gz
dobase=${fname%_1.clean.fq.gz}sample=${fname%%_*}bwa mem -t 20 -M -R "@RG\tID:${sample}\tSM:${sample}\tPL:illumina\tLB:${sample}" ../Ref-T2T/T2T_genome.fa "$fname" "${base}_2.clean.fq.gz" |samtools view -b -@ 5 -|samtools sort -m 5g -@ 5 - > "../1.sortbam/${base}_sort.bam"
done

##代码解释：
../Ref-T2T/T2T_genome.fa表示参考基因组文件的路径。(按具体情况修改)"$fname"表示输入的clean.fq.gz文件的路径。"${base}_2.clean.fq.gz"表示对应的_2.clean.fq.gz文件的路径。|samtools view -b -@ 5 -表示将比对结果从SAM格式转换为BAM格式。|samtools sort -m 5g -@ 5 -表示对BAM文件进行排序。"../1.sortbam/${base}_sort.bam"表示排序后的bam文件的输出路径。（按具体情况修改）

7. 按个体文件编号统计sort.bam文件

depth.txt是测序深度，coverage.txt是覆盖度，coverage_4X.txt是4x覆盖度

这段代码主要用于统计和分析 bam 文件中的比对结果和覆盖度信息

工作目录：在前一步创建的1.sortbam目录下运行

#运行此脚本
for filename in *_sort.bam
do
samtools flagstat $filename|paste -sd "/t" - >> stat_res.txt
samtools depth $filename|awk '{sum+=$3}END{print sum/NR}' >> depth.txt
samtools depth $filename|awk '{if($3>0)a+=1}END{print a/NR}' >> coverage.txt
samtools depth $filename|awk '{if($3>4)a+=1}END{print a/NR}' >> coverage_4X.txt
done

stat_res.txt文件解读：

上述结果是按照bam文件在文件夹中的排序得到的：
（1）获得所有bam文件的id

ls *bam>bam.id.txt

（2）合并个体编号的stat的结果

paste bam.id.txt stat_res.txt >248stat_res.txt #用excel提取总total 和mapped的 reads数，以及mapped的比例，

（3）利用vi 给depth.txt、coverage.txt 和coverage_4X.txt加抬头，并合并

paste depth.txt coverage.txt coverage_4X.txt >248.txt

（4）去掉excel处理后的文件中的^M，然后再和上述结果合并

sed -e 's/.$//' 248stat_res.txt >248stat_res.new.txt    #将每行文本的最后一个字符删除
paste 248stat_res.new.txt 248.txt >248tongji.txt

补充：用代码提取总total 和mapped的 reads数，以及mapped的比例，（以下代码皆为本人使用，至于是否具有普适性，请自行修改和判断）

#1.
import pandas as pd
#按照文件中第二列内容中的“/”号分割为多列（N/A不分割）
df = pd.read_csv(r"C:\Users\Desktop\248stat_res.txt", delimiter="\t", header=None)
# 提取第一列数据
column1 = df.iloc[:, 0]
# 将第二列按 "/" 分割成多列
df_split = df.iloc[:, 1].str.split('/', expand=True)
# 创建一个空的列表用于存储处理后的列
processed_columns = []
# 合并相邻的首尾为 "N" 或 "A" 的列
i = 0
while i < len(df_split.columns):if i < len(df_split.columns) - 1 and (df_split[i].str.endswith('N') | df_split[i].str.endswith('A')).all() and (df_split[i + 1].str.startswith('N') | df_split[i + 1].str.startswith('A')).all():processed_columns.append(df_split[i] + '/' + df_split[i + 1])i += 2else:processed_columns.append(df_split[i])i += 1
# 将处理后的列合并成新的DataFrame
df_processed = pd.concat([column1, pd.concat(processed_columns, axis=1)], axis=1)
# 打印处理后的数据
# print(df_processed)
# df_l = df_processed.iloc[1,:]
# print(len(df_l))
# for i in df_l:
#     print(i)
# 将处理后的DataFrame保存为txt文件
df_processed.to_csv(r"C:\Users\Desktop\processed_248stat_res.txt", sep='\t', index=False, header=False)#2.#用excel提取总total 和mapped的 reads数，以及mapped的比例，#即形成四列（Sample、total reads、mapped reads、mapping rate)# 读取txt文件
df = pd.read_csv(r"C:\Users\Desktop\processed_248stat_res.txt", delimiter="\t", header=None)
Sample = df.iloc[:, 0]
print(Sample)
# 提取第二列和第五列数据
column2 = df.iloc[:, 1]
column5 = df.iloc[:, 4]
# # 打印提取的数据
# print("第二列数据:")
# print(column2)
#
# print("\n第五列数据:")
# print(column5)# 使用正则表达式提取数字
total_reads = column2.str.extract(r'(\d+)')[0]
# print(total_reads)# 提取第五列数据中的数字和百分比
extracted_data = column5.str.extract(r'duplicatest(\d+) \+ \d+ mapped \((\d+\.\d+)% :')# 添加列标签
extracted_data.columns = ['mapped reads', 'mapping rate']# 添加百分号到 Mapping Rate 列
extracted_data['mapping rate'] = extracted_data['mapping rate'] + '%'extracted_data['Sample'] = Sample
extracted_data['total reads'] = total_reads
extracted_data =extracted_data[['Sample','total reads','mapped reads','mapping rate']]
# 打印提取的数据
print(extracted_data)
# 将处理后的extracted_data保存为txt文件
extracted_data.to_csv(r"C:\Users\Desktop\finally_248stat_res.txt", sep='\t', index=False)#3.
# 按照另一个表格文件的样本名相应的将本txt文件中对应的内容填写（将填写的表格内容复制即可）
df = pd.read_csv(r"C:\Users\Desktop\248tongji.txt", delimiter="\t")
# print(df)
Sample = df.iloc[:, 0]
split_Sample = Sample.str.split("_")
first_elements = [item[0] for item in split_Sample]
print(first_elements)
df['Sample'] = first_elements
print(df)
# print(new_Sample)
#读取xlsx文件中的一个表（248比对到T2T统计）
# 读取xlsx文件
file_path = r"C:\Users\Desktop\248重测序样本20240706.xlsx"  
df_xlsx = pd.read_excel(file_path, sheet_name='248比对到T2T统计')
# 打印读取的表格内容
print(df_xlsx)
# 根据第一列Sample列的样本数据查询并填充信息
df_merged = df.merge(df_xlsx, on="Sample", how="right")
print(df_merged)
df_merged.to_csv(r"C:\Users\Desktop\248比对到T2T统计.txt", sep='\t', index=False)

8. 简单过滤，去除unmapped和mate unmapped的reads生成sort_flt.bam文件

先在上一级目录创建2.sort_flt.bam文件夹
代码工作目录：在1.sortbam目录下运行

#运行guolv.sh脚本：
for fname in *_sort.bam
do
sample=${fname%_sort.bam}
samtools view -q 20 -f 0x0002 -F 0X0004 -F 0X0008 -b $fname
>../2.sort_flt.bam/${sample}_sort_flt.bam
done

9.标记重复生成srt_flt_dup.bam文件

先在上一级目录创建3.srt_flt_dup.bam文件夹

代码工作目录：在上面创建的2.sort_flt.bam目录下运行

#运行bioajichongfu.sh脚本
for fname in *_sort_flt.bam
do
sample=${fname%_sort_flt.bam}
picard MarkDuplicates I=$fname O=../3.srt_flt_dup.bam/${sample}_srt_flt_dup.bam
M=../3.srt_flt_dup.bam/${sample}.duplicates.log
done

10. 使用GATK的HaplotypeCaller进行SNP/INDEL分析

先在上一级目录创建4.srt_flt_dup.g.vcf文件夹

工作路径：上一步创建的3.srt_flt_dup.bam文件夹

#先对文件夹中的bam文件建索引
for fname in *_srt_flt_dup.bam
do
samtools index $fname
done

#再使用gatk HaplotypeCaller进行SNP/INDEL分析
for fname in *_srt_flt_dup.bam
do
sample=${fname%_srt_flt_dup.bam}
gatk HaplotypeCaller -R ../Ref-T2T/T2T_genome.fa -I $fname -ERC GVCF -O
../4.srt_flt_dup.g.vcf/${sample}_srt_flt_dup.g.vcf.gz
done

11.使用gatk的CombineGVCFs命令将前面鉴定的坛紫菜样品的.g.vcf文件整合到一起

这里for循环不好用，需要用-V 把每个文件依次输入。最终获得248toTZC.g.vcf.gz文件
先在上一级目录创建5.haitanensis文件夹

工作路径：上一步创建的4.srt_flt_dup.g.vcf文件夹

可以直接用-V 来一个一个添加待整合的vcf文件（但是太慢了，可以考虑使用脚本来循环处理）

为了方便，可以在当前目录下创建一个combine.sh脚本，用于合并vcf文件（注意：如果是之前第4步已经挑选了所需要的坛紫菜样品，则直接获取当前目录中所有匹配的 VCF 文件，也就是用以下的代码）

#!/bin/bash# 指定参考基因组文件
REF_GENOME="../Ref-T2T/T2T_genome.fa"# 指定输出文件
OUTPUT="../5.haitanensis/248toTZC.g.vcf.gz"# 获取当前目录中所有匹配的 VCF 文件
VCF_FILES=$(ls *_srt_flt_dup.g.vcf.gz)# 初始化 -V 参数
V_ARGS=""# 遍历 VCF 文件列表
for VCF_FILE in $VCF_FILES; do# 获取文件名（不包括扩展名）FILE_BASE=$(basename $VCF_FILE _srt_flt_dup.g.vcf.gz)# 为每个文件添加 -V 参数V_ARGS="$V_ARGS -V $FILE_BASE\_srt_flt_dup.g.vcf.gz"
done# 使用 CombineGVCFs 合并 VCF 文件
gatk CombineGVCFs -R $REF_GENOME $V_ARGS -O $OUTPUT

如果之前第4步没有立即挑选出所需要的坛紫菜样品，而是只是用一个表格记录下来；那么要从当前目录中获取坛紫菜样品的vcf文件。使用以下代码：

#!/bin/bash
# 指定参考基因组文件
REF_GENOME="../Ref-T2T/T2T_genome.fa"
# 指定输出文件
OUTPUT="../5.haitanensis/149toTZC.g.vcf.gz"
# 获取当前目录中所有匹配的 VCF 文件
VCF_FILES=$(cat tang_fixed.txt | xargs -n 1 ls)
# 初始化 -V 参数
V_ARGS=""# 遍历 VCF 文件列表
for VCF_FILE in $VCF_FILES; doif [ -f "$VCF_FILE" ]; thenecho "文件存在: $VCF_FILE"# 获取文件名（不包括扩展名）FILE_BASE=$(basename $VCF_FILE _srt_flt_dup.g.vcf.gz)# 为每个文件添加 -V 参数V_ARGS="$V_ARGS -V $FILE_BASE\_srt_flt_dup.g.vcf.gz"elseecho "文件不存在: $VCF_FILE"fi
done
# 使用 CombineGVCFs 合并 VCF 文件
gatk CombineGVCFs -R $REF_GENOME $V_ARGS -O $OUTPUT

补充：（tang_fixed.txt文件形成）
在这之前需要将cp构树鉴定结果为“坛紫菜”的挑选出来，获取对应重测序分析编号，以便获取对应的vcf文件名称。

#使用代码形成为坛紫菜的vcf文件名称，保存到tang.txt文件中
import pandas as pd
file_path = r"C:\Users\Desktop\tang.xlsx"
df_xlsx= pd.read_excel(file_path)
# print(df_xlsx)
print(df_xlsx.columns)
# print(df_xlsx.index)
# # 筛选数据
filtered_data = df_xlsx[df_xlsx["cp构树鉴定结果"] == "坛紫菜"]
# 重置索引
filtered_data.reset_index(drop=True, inplace=True)
# 输出结果
print(filtered_data)
# 获取第一列数据
first_column = filtered_data.iloc[:, 0]
# 将第一列数据写入文本文件
with open("tang.txt", "w") as f:for item in first_column:f.write(str(item) + "_srt_flt_dup.g.vcf.gz\n")

如果文件名之间有分隔符，可能会导致问题。

#删除 tang.txt 文件中的所有 \r 字符，然后将结果保存到名为tang_fixed.txt 的新文件
tr -d '\r' < tang.txt > tang_fixed.txt

12. 获得原始vcf文件

工作路径：上一步创建的5.haitanensis文件夹
运行此脚本

echo 'start_get_rawVCF'
gatk GenotypeGVCFs -R ../Ref-T2T/T2T_genome.fa -O
149toTZC.variant.raw.vcf.gz -V 149toTZC.g.vcf.gz -all-sites 1>
log_GenotypeGVCFs_all_sites.txt 2>&1
echo 'done_get_rawVCF'

13.保留双等位基因型

工作路径：5.haitanensis文件夹

vcftools --gzvcf 149toTZC.variant.raw.vcf.gz --min-alleles 2 --max-alleles 2 --
recode --recode-INFO-all --out 149toTZC.2allele.vcf

14.最终得到vcf文件

查看一下得到的文件

less -N 149toTZC.2allele.vcf.recode.vcf

sed -n '40,70p' 149toTZC.2allele.vcf.recode.vcf | cut -f 1-6,10-14

---

总结

本文主要了对重测序数据进行处理得到vcf文件的流程，从得到rawdata，用fastp对数据进行质控和过滤，使用getorganelle软件对叶绿体基因组进行组装、通过构树鉴定所需个体；
再到后期的给参考基因组建索引、clean reads比对到参考基因组生成bam文件、对bam文件进行一系列处理、使用gatk软件进行SNP/INDEL分析、文件整合、获得原始vcf文件。

总之，得到了vcf文件后，是用于后续的继续分析。

2024/7/21