群体重测序数据分析笔记

在生物信息学中，群体重测序数据的挖掘和分析对于理解生物的进化、自然选择以及功能基因的定位等研究具有重要的意义。今天分享的笔记内容是群体遗传学相关的知识点，下面将一步步介绍整个重测序分析的流程和方法。

分析常用流程和方法简述

常用的重测序分析流程一般包含以下步骤：

1. 质控和数据准备：

这一步包括对原始测序数据进行质量评估和控制，以及数据格式的转换，常用流程包括测序数据的质控(QC)、比对、标记重复、排序、建立索引和变异检测。常用的软件工具包括FastQC、BWA、SAMtools、Picard和GATK等。

首先，我们需要进行质量控制，确保我们的数据是准确可靠的。常用的质量控制工具有FastQC和Trimmomatic。以下是用Trimmomatic进行质量控制的例子：

java -jar trimmomatic.jar PE -phred33 \
input_forward.fq.gz input_reverse.fq.gz \
output_forward_paired.fq.gz output_forward_unpaired.fq.gz \
output_reverse_paired.fq.gz output_reverse_unpaired.fq.gz \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 \
SLIDINGWINDOW:4:15 MINLEN:36

2. 序列比对：

将处理好的数据比对到参考基因组上，得到比对结果,比对可以使用BWA这个工具。以下是比对命令的示例：

# 建立参考基因组索引
bwa index ref.fa
# 比对
bwa mem ref.fa read1.fq read2.fq > aln.sam

3. 变异检测：

根据比对结果，检测和分析基因序列中的变异,变异检测通常使用SAMtools和BCFtools。以下是如何使用这些工具进行变异检测的示例：

# 转换格式
samtools view -S -b aln.sam > aln.bam
# 排序
samtools sort aln.bam -o aln_sorted.bam
# 检测变异
samtools mpileup -uf ref.fa aln_sorted.bam | bcftools call -cv - > var.raw.vcf

4. 变异注释和分析：

对检测到的变异进行注释和深度分析，包括群体结构分析、选择性消除分析、全基因组关联分析等。

以下为流程简单示意

# 质量控制
fastqc raw_data.fq
# 比对
bwa mem reference.fa raw_data.fq > aln.sam
# 标记重复
picard MarkDuplicates I=aln.sam O=marked.bam M=metrics.txt
# 排序
samtools sort -O bam -o sorted.bam -T temp.prefix marked.bam
# 建立索引
samtools index sorted.bam
# 变异检测
gatk HaplotypeCaller -R reference.fa -I sorted.bam -O variants.vcf

上游分析变异检测方法代码

在进行测序数据序列文件上游变异检测时，我们通常使用如GATK工具。以下是一段使用GATK进行SNP和Indel检测的代码：

# SNP检测
gatk --java-options "-Xmx4g" \
HaplotypeCaller  -R reference.fa \
-I sorted.bam \
-O raw_snps.vcf
# Indel检测
gatk --java-options "-Xmx4g" \
HaplotypeCaller  -R reference.fa \
-I sorted.bam --ploidy 2 \
--genotyping-mode DISCOVERY \
-stand_emit_conf 10 \
-stand_call_conf 30 \
-O raw_indels.vcf

群体结构分析方法

对于群体结构的分析，我们通常会使用程序如PLINK和ADMIXTURE。PLINK可以用于生成适合ADMIXTURE分析的数据，而ADMIXTURE可以进行群体结构分析。

# PLINK生成适用于ADMIXTURE的数据
plink --file input_data \
      --make-bed --out plink_output
# ADMIXTURE进行群体结构分析
admixture --cv plink_output.bed K

变异位点选择性消除分析

在基因组变异位点选择性消除分析中，我们可以使用vcftools和R等工具来评估基因位点多态性和选择分化差异。以下是具体的操作：

# 使用vcftools计算窗口内的多态性
vcftools --vcf var.raw.vcf --window-pi 50000

# 在R中进行选择分化差异判断
# 以下是示例代码，具体代码需要根据实际情况编写
library(ape)
data <- read.table("fst.txt", header = T)
fst <- data$V4
sig_level <- qnorm(1 - 0.05 / 2) / sqrt(2)
outliers <- which(fst > sig_level)

GWAS全基因组关联分析

GWAS全基因组关联分析我们可以使用GAPIT包进行分析。以下是用R语言进行GWAS全基因组关联分析的示例代码：

library(GAPIT)
genotype_file <- "genotype.hmp.txt"
phenotype_file <- "phenotype.txt"
GAPIT_data <- GAPIT.Data(fileHapmap = genotype_file,
                         filePhenotype = phenotype_file)
GAPIT_GLM <- GAPIT.GLM(GAPIT_data)
GAPIT_MLM <- GAPIT.MLM(GAPIT_data)

PCA分析与进化树绘制

我们可以使用plink和gcta进行PCA分析，并使用ggtree进行群体结构进化树的绘制。以下是具体的操作：

# 使用plink和gcta进行PCA分析
plink --bfile plink --make-grm-bin --out plink
gcta --grm-bin plink --pca 10 --out plink

# 在R中使用ggtree绘制进化树
# 以下是示例代码，具体代码需要根据实际情况编写
library(ggtree)
tree <- read.tree("tree.nwk")
ggtree(tree) + geom_tiplab()

从vcf文件中挖掘显著变异位点的频率变化信息

我们可以使用vcftools工具从vcf文件中挖掘显著变异位点的频率变化信息，以下是具体的操作：

vcftools --vcf var.raw.vcf --freq --out freq

以上就是群体重测序数据的挖掘与分析的整个过程，每一步都需要我们仔细和认真的处理。

总结

以下是一个简单的bash脚本，实现对多个样品的批量重测序分析,这个脚本使用了一个循环来处理每一个样品。请注意这只是一个示例脚本，您可能需要根据实际情况对其进行修改或优化。

#!/bin/bash

# 路径参数
REF_GENOME_PATH=/path/to/your/reference/genome
SAMPLE_LIST=/path/to/your/sample/list
SAMPLE_PATH=/path/to/your/sample/data
WORKING_DIR=/path/to/your/working/directory

# 创建工作目录
mkdir -p ${WORKING_DIR}

# 循环处理每一个样品
while read SAMPLE; do
    echo "Processing sample ${SAMPLE}..."

    # 1. 质控和数据准备
    java -jar trimmomatic.jar PE -phred33 \
        ${SAMPLE_PATH}/${SAMPLE}_1.fq.gz \
        ${SAMPLE_PATH}/${SAMPLE}_2.fq.gz \
        ${WORKING_DIR}/${SAMPLE}_1_paired.fq.gz \
        ${WORKING_DIR}/${SAMPLE}_1_unpaired.fq.gz \
        ${WORKING_DIR}/${SAMPLE}_2_paired.fq.gz \
        ${WORKING_DIR}/${SAMPLE}_2_unpaired.fq.gz \
        ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

    # 2. 序列比对
    bwa mem ${REF_GENOME_PATH} \
        ${WORKING_DIR}/${SAMPLE}_1_paired.fq.gz \
        ${WORKING_DIR}/${SAMPLE}_2_paired.fq.gz > ${WORKING_DIR}/${SAMPLE}.sam

    # 3. 变异检测
    samtools view -S -b ${WORKING_DIR}/${SAMPLE}.sam > ${WORKING_DIR}/${SAMPLE}.bam
    samtools sort ${WORKING_DIR}/${SAMPLE}.bam -o ${WORKING_DIR}/${SAMPLE}_sorted.bam
    samtools mpileup -uf ${REF_GENOME_PATH} ${WORKING_DIR}/${SAMPLE}_sorted.bam | bcftools call -cv - > ${WORKING_DIR}/${SAMPLE}_var.raw.vcf

    # 4. 提取VCF文件中的频率信息
    vcftools --vcf ${WORKING_DIR}/${SAMPLE}_var.raw.vcf --freq --out ${WORKING_DIR}/${SAMPLE}_freq

    echo "Finished processing sample ${SAMPLE}."

done < ${SAMPLE_LIST}

在这个脚本中，我们假设有一个文本文件 ${SAMPLE_LIST} 包含所有待处理的样品名称，每个样品对应一对fastq.gz文件，文件名分别为 ${SAMPLE}_1.fq.gz 和 ${SAMPLE}_2.fq.gz。所有这些文件都存储在 ${SAMPLE_PATH} 路径下。处理过程中的中间文件和结果文件都会存储在 ${WORKING_DIR} 路径下。

使用这个脚本之前，你需要修改上面的四个路径参数，使它们指向实际的路径。另外，这个脚本只执行了一部分的分析步骤，如果你需要执行更多的步骤，你可以在脚本中添加相应的命令。