当前位置：首页 > article >正文

BSR脚本

article 2025/1/31 12:37:19

1. 参数输入
   └── 输入 FASTQ 文件（突变体组和野生型组）

2. 比对及生成 BAM 文件
   ├── HISAT2 或 Tophat2 构建参考基因组索引
   ├── 使用 HISAT2 或 Tophat2 进行比对（生成 SAM 文件）
   └── 使用 Samtools 转换 SAM 文件为 BAM 文件

3. 读取分组并对 BAM 文件排序
   ├── 使用 Picard 添加读取组信息
   ├── 使用 Samtools 排序 BAM 文件
   ├── 使用 Picard ReorderSam 调整染色体顺序
   ├── 使用 Samtools 为排序后的 BAM 文件生成索引
   └── 重复处理野生型组的 BAM 文件

4. SNP 检测
   ├── 使用 GATK 的 UnifiedGenotyper 进行 SNP 和 InDel 检测
   └── 生成原始的 SNP 和 InDel VCF 文件

5. SNP 过滤
   ├── 使用 GATK 的 SelectVariants 进行硬过滤
   ├── 应用质量过滤标准（QD、FS、MQRankSum 等）
   └── 输出筛选后的 VCF 文件

6. 压缩与索引 VCF 文件
   ├── 使用 bgzip 压缩 VCF 文件
   └── 使用 tabix 创建 VCF 文件的索引

7. SNP 绘图
   ├── 通过 BSR 分析生成 SNP 分布数据
   ├── 使用 R 脚本绘制 SNP 分布图
   └── 输出结果图表（PDF 格式）

第1部分：参数输入（可更改）

mut_read1=$1 # 变异体组的fastq文件1（gzip压缩）
mut_read2=$2 # 变异体组的fastq文件2（gzip压缩）
mut_line_name=$3 # 变异体组名称

wt_read1=$4 # 野生型组的fastq文件1（gzip压缩）
wt_read2=$5 # 野生型组的fastq文件2（gzip压缩）
wt_line_name=$6 # 野生型组名称

第2部分：比对及重命名TopHat/hisat2生成的BAM文件

参考基因组路径

# reference=/public/home/04034/Maize_reference/B73-ensembl-V4-50
reference=~/yiyaoran/workspace/06.BSRseq/guo_BSR_pipline/1516/ref/genome.fasta

加载GCC模块

module load compiler/gcc/7.3.1

方法1.使用 HISAT2

1.1构建参考基因组索引

hisat2-build <reference.fasta> <output_index_prefix>

运行成功后，HISAT2 将生成 6 个索引文件，文件名格式为：

genome_index.1.ht2
genome_index.2.ht2
genome_index.3.ht2
genome_index.4.ht2
genome_index.5.ht2
genome_index.6.ht2

1.2使用Hisat2对变异体组进行比对，并将输出的SAM文件转换为BAM文件

# 使用 HISAT2 对变异样本的成对读段进行比对
hisat2 -x $reference \                 # 指定参考基因组索引文件的前缀（通常为 .ht2 文件）
    -1 $mut_read1 -2 $mut_read2 \      # 指定用于比对的成对 FASTQ 文件（$mut_read1 和 $mut_read2 分别为 read1 和 read2）
    -S $mut_line_name".sam" \          # 指定输出的 SAM 文件名称，以样本名称为前缀 ($mut_line_name)
    --new-summary 2>$mut_line_name".hisat2.log"  # 将新的比对摘要日志信息重定向输出到日志文件中

# 使用 samtools 将 SAM 文件转换为 BAM 格式
samtools view -bS $mut_line_name".sam" > $mut_line_name".bam"  # 将 SAM 文件转换为 BAM 文件，输出命名为样本前缀 ($mut_line_name)



# 对野生型组进行相同的比对和转换
hisat2 -x $reference -1 $wt_read1 -2 $wt_read2 -S $wt_line_name".sam" --new-summary 2>$wt_line_name".hisat2.log"
samtools view -bS $wt_line_name".sam" > $wt_line_name".bam"

echo "###########################"
echo "# 比对及重命名完成 #"
echo "###########################"

.bam文件

QNAME       FLAG    RNAME      POS    MAPQ    CIGAR    RNEXT    PNEXT    TLEN     SEQ                                                                    QUAL                                                                       AS:i   ZS:i  XN:i  XM:i  XO:i  XG:i  NM:i  MD:Z       YS:i  YT:Z  NH:i
E200016962L1C001R00100012647  83     scaf_258 49446   1       150M     =        14179   -35417   GGCCGCCTCTCCAGACTACAATTCGGGCGGCGAGGCCGCCCGATTCTCAAGCTGGGCTGCTCCCGGTTCGCTCGCCGTTACTAGGGGAATCCTCGTAAGTTTCTTCTCCTCCGCTTATTTATATGCTTAAACTCAGCGGGTAGTCCCGAC  GC"GFGGAFFFFEFF<F?GFAFF8GF4FFGFFFFFF<FEFFFFFFFFFGFFFFEFFF?FFFFFFFFFFFFFFFFFFFFFFFFFGGFFFFFFFFFFFFDFFFEFFCFFFFFFEFFFFFFFFFFFFFGFFFFFFEFFFFFFFFFFFF>FFFF  AS:i:-2 ZS:i:-2 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:2A147 YS:i:0 YT:Z:CP NH:i:4
E200016962L1C001R00100012647  163    scaf_258 14179   1       150M     =        49446   35417    GCAATTTCAAGCACTCTTTGACTCTCTTTTCAAAGTCCTTTTCATCTTTCCCTCGCGGTACTTGTTCGCTATCGGTCTCTCGCCTGTATTTAGCCTTGGACGGAGTTTACCGCCCGATTTGGGCTGCATTCCCAAACAACCCGACTCGTT  FFFFFGAFFFFFFFGFGFFFEFGFGFGGFGFFFFFGFFGDFFFFFFGFGFFFGFEEGFGFFGGFFFFFFGFGFFFFFGFFFFFFGGFEFGFFFFFFGGFCFFF;FFGDGFEFFFFFAGFGFFDFGFFFGGFFGFFFFFFFGGFFFGFFGG  AS:i:0  ZS:i:0  XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:150   YS:i:-2 YT:Z:CP NH:i:4
1. QNAME:    Query template/pair NAME（查询模板或配对的名称）
2. FLAG:     Bitwise FLAG（标志，用于表示比对结果的属性）
3. RNAME:    Reference sequence NAME（参考序列的名称，即比对到的染色体或 scaffold 名称）
4. POS:      1-based leftmost mapping POSition（比对的起始位置，以1为基准）
5. MAPQ:     Mapping Quality（比对质量得分）
6. CIGAR:    CIGAR string（表示比对过程中的匹配、插入、缺失等情况）
7. RNEXT:    Reference name of the mate/next read（配对 read 或下一个 read 的参考序列名称）
8. PNEXT:    Position of the mate/next read（配对 read 或下一个 read 的比对位置）
9. TLEN:     Observed Template LENgth（观测到的模板长度）
10. SEQ:     Segment SEQuence（序列，比对的 read 序列）
11. QUAL:    ASCII of Phred-scaled base QUALity+33（碱基质量得分）
12. TAGS:    Optional fields in format TAG:TYPE:VALUE（可选字段，包含额外的比对信息）

示例解释：

E200016962L1C001R00100012647	83	scaf_258	49446	1	150M	=	14179	-35417	GGCCGCCTCTCCAGACTACAATTCGGGCGGCGAGGCCGCCCGATTCTCAAGCTGGGCTGCTCCCGGTTCGCTCGCCGTTACTAGGGGAATCCTCGTAAGTTTCTTCTCCTCCGCTTATTTATATGCTTAAACTCAGCGGGTAGTCCCGAC	GC"GFGGAFFFFEFF<F?GFAFF8GF4FFGFFFFFF<FEFFFFFFFFFGFFFFEFFF?FFFFFFFFFFFFFFFFFFFFFFFFFGGFFFFFFFFFFFFDFFFEFFCFFFFFFEFFFFFFFFFFFFFGFFFFFFEFFFFFFFFFFFF>FFFF	AS:i:-2	ZS:i:-2	XN:i:0	XM:i:1	XO:i:0	XG:i:0	NM:i:1	MD:Z:2A147	YS:i:0	YT:Z:CP	NH:i:4

列名	示例数据	解释
QNAME	E200016962L1C001R00100012647	查询序列的名称（read 的名称或标识符）
FLAG	83	比对标志，表明 read 的比对情况及状态
RNAME	scaf_258	比对到的参考序列名称（染色体、contig 或 scaffold）
POS	49446	比对的起始位置
MAPQ	1	比对质量得分，通常范围是 0-60，数值越高表示比对越可靠
CIGAR	150M	CIGAR 字符串，表示 150 个连续匹配（150 bp）
RNEXT	=	配对 read 的参考序列名称（`=` 表示与 RNAME 相同）
PNEXT	14179	配对 read 的比对起始位置
TLEN	-35417	观测到的模板长度（负值表示两个配对 read 之间的方向）
SEQ	GGCCGCCTCTCC…	比对的序列
QUAL	GC"GFGGAFFFFE…	碱基质量得分（Phred 编码，表示每个碱基的测序质量）
TAGS	AS:i:-2 ZS:i:-2 XN:i:0 XM:i:1 …	可选字段，包括比对的得分、错误信息和比对状态等

TAGS 字段解释：

AS:i: 比对得分（alignment score）
ZS:i: 二次比对得分（second alignment score）
XN:i: 非匹配的基因组位置数量
XM:i: 碱基错配数（number of mismatches）
XO:i: 碱基插入/缺失的操作数（number of gap operations）
XG:i: 碱基插入/缺失的长度（gap length）
NM:i: 碱基错配、插入或缺失的总数（edit distance）
MD:Z: MD 字符串，用于指示比对中的错配位置及替换的碱基
YS:i: 配对 read 的比对得分
YT:Z: 比对的类型
NH:i: 比对到多个位置的次数（number of hits）

方法2使用Tophat2进行

tophat2 -p 4 \  # 使用 4 个线程进行并行计算，加快比对速度
    -g 1 \  # 每个 read 只允许比对到基因组中的一个位置，减少多拷贝基因的干扰
    --mate-inner-dist 100 \  # 配对测序中两个 reads 之间的平均插入片段距离设置为 100 碱基
    --mate-std-dev 75 \  # 配对 reads 插入片段长度的标准差设置为 75 碱基
    --microexon-search \  # 启用微外显子搜索功能，以便检测较短的外显子
    -o $mut_line_name $reference  # 设置输出目录名称为 $mut_line_name，使用的参考基因组为 $reference

第3部分：读取分组并对BAM文件排序

samtools sort：用于对 BAM 文件中的 reads 进行位置排序，以便下游分析更高效。
ReorderSam：用于调整染色体的排列顺序，以确保它与参考基因组完全一致。

3.1使用 Picard 工具为 BAM 文件添加或替换读取组信息

java -jar ~/yiyaoran/guoyuhang/software/picard-tools-1.119/AddOrReplaceReadGroups.jar \
    I=$mut_line_name.bam \  # 输入 BAM 文件
    O=$mut_line_name"_GR.bam" \  # 输出 BAM 文件，添加或替换读取组后的文件名
    RGID=$mut_line_name \  # 读取组 ID（Read Group ID），一般使用样本名称
    RGLB=A \  # 文库 ID（Library），此处设置为 A
    RGPL=illumina \  # 测序平台（Platform），例如 Illumina
    RGPU=A \  # 测序单元（Platform Unit），此处设置为 A
    RGSM=$mut_line_name  # 样本名称（Sample），与样本名称一致

生成文件如下：_GR.bam内容如下:
与.bam区别：6463-1z_GR.bam中，每条记录都包含了一个 RG:Z:6463-1z 标签，这个标签表示“Read Group”（读取组）。这个标签用于标识不同样本、文库、或测序批次的信息。

samtools view 6463-1z_GR.bam | less -S | head -n 2
E200016962L1C001R00100012647	83	scaf_258	49446	1	150M	=	14179	-35417	GGCCGCCTCTCCAGACTACAATTCGGGCGGCGAGGCCGCCCGATTCTCAAGCTGGGCTGCTCCCGGTTCGCTCGCCGTTACTAGGGGAATCCTCGTAAGTTTCTTCTCCTCCGCTTATTTATATGCTTAAACTCAGCGGGTAGTCCCGAC	GC"GFGGAFFFFEFF<F?GFAFF8GF4FFGF