跳转至

基因组重测序

在生物信息学中，变异检测（Variant Calling） 是一个关键步骤，旨在识别基因组中的遗传变异（如单核苷酸多态性（SNP）和插入缺失（InDel））。这些变异可以提供关于疾病、进化和个体特征的宝贵信息。

一 准备样品信息

在分析前，需要整理和准备样品的基础信息，包括以下内容：

  • 样品名称：每个样品的唯一标识。
  • FASTQ 文件路径：对应的测序文件路径。
  • 测序平台：例如 Illumina 或其他平台。
  • 是否为 Clean Data：标注数据是否已经过预处理。

这些信息将为后续分析流程的调用和自动化处理提供便利。

Unit Sheet 文件

为了统一和规范样品的管理，需要为每个样品创建一个或多个测序单元（unit），并将其信息记录在 units.txt 文件中。文件格式如下：

列名 描述
sample 样品名称，表示生物样品的唯一标识
fq1, fq2 双端测序数据的 FASTQ 文件路径
其他列（可选） 额外的信息，例如测序平台、数据质量等

文件格式要求

在 units.txt 文件中，每一行表示一个样品的测序单元，必须包含以下内容：

  1. 样品名称：通过 sample 列与生物样品关联。
  2. FASTQ 文件路径
  3. fq1 和 fq2 列分别指向双端测序的两个 FASTQ 文件。
  4. 如果某些值缺失，可通过空列指定。

示例文件内容

sample  fq1                      fq2
S1      /path/to/S1_R1.fq.gz    /path/to/S1_R2.fq.gz
S2      /path/to/S2_R1.fq.gz    /path/to/S2_R2.fq.gz
S3      /path/to/S3_R1.fq.gz

通过这样的样品信息表格，可以轻松管理样品并与分析流程对接，确保高效、准确的处理。

生成 units 表格

echo -e "sample\tfq1\tfq2" > units.txt
ls *.gz |awk -F'_' '{print $1}' |sort -u |while read sample ;
do 
  fq1=$(find $(pwd) -name "${sample}*1*.gz")
  fq2=$(find $(pwd) -name "${sample}*2*.gz")
  echo -e "$sample\t$fq1\t$fq2";
done >> units.txt

可以通过以下方式方便的获取测序文件信息，进行相关操作。

sed '1d' units.txt |while read sample fq1 fq2;
do
  # 在此处对每个样本进行分析操作
  echo "Processing sample: $sample"
  echo "Fastq file 1: $fq1"
  echo "Fastq file 2: $fq2"
done

二 数据质控

为保证下游分析输入数据的可靠性，需要对下机原始测序数据进行质控。目前，测序公司通常会提供clean data，一般情况无需额外质控，如需质控一般会使用以下工具。

  1. Fastqc
  2. Multiqc
  3. Fastp

三 Variant Calling

下面是一些常用的call变异流程

bwa + GATK

# LSF 集群
tail -n +2 data/units.txt |while read i fq1 fq2;
do 
    bsub -J ${i}_gatk_1 -n 5 -q q2680v2 -R "span[hosts=1]" -o $i.callVCF.out -e $i.callVCF.err \
        sh 1.gatk_calling.sh $i $fq1 $fq2
    sleep 3
done

# 加载所需软件
module load Singularity/3.7.3
module load SAMtools/1.9
module load BWA/0.7.17
module load GATK/4.5.0.0
module load picard/2.23.9
module load R/4.0.0
sif=~/singularity.sif/reseq.sif
# 样品相关变量获取
sample=$1
i=$sample
fq1=$2
fq2=$3

#设置参考基因相关文件变量，方便后续使用
refdir=~/fast/reseq/ref
REF=$refdir/Arabidopsis_thaliana_dna.Chr4.fa
BWA_INDEX=$refdir/index/bwa-mem2/Arabidopsis_thaliana_dna.Chr4.fa
GFF=$refdir/Arabidopsis_thaliana.Chr4.gff3
GTF=$refdir/Arabidopsis_thaliana.Chr4.gtf
INDEX_FAI=$refdir/Arabidopsis_thaliana_dna.Chr4.fa.fai
PICARD_DICT=$refdir/Arabidopsis_thaliana_dna.Chr4.dict

# ******************************************
# 0. Setup
# ******************************************
# 需要使用的核心数
nt=5
# 比对信息
group="G"
platform="ILLUMINA"
mq=30

# 工作目录
workdir=~/fast/reseq/result/
[ ! -d $workdir ] && mkdir -p $workdir
tmpdir=$workdir/temp
[ ! -d $tmpdir ] && mkdir -p $tmpdir
cd $workdir
# 输出文件
rawbam=$i.bam
sortedbam=$i.sorted.bam
depbam=$i.deduped.bam
outvcf=${i}_raw.gvcf
exec > $workdir/$i.callVCF.log 2>&1 # call vcf的日志文件

# ******************************************
# 1. 利用 BWA-MEM2 进行比对并排序
# ******************************************
echo -e "\n\nStep 1: BWA-MEM2 mapping and sorting"
(bwa-mem2 mem -t ${nt} -M -R "@RG\tID:${group}\tLB:${i}\tPL:${platform}\tSM:${i}" \
    $BWA_INDEX  ${fq1} ${fq2} || echo -n 'error') \
    | samtools sort -@ ${nt} > $rawbam && samtools index -@ $nt $rawbam 
samtools view -hbS  -q $mq -o $sortedbam $rawbam && \
samtools index -@ $nt $sortedbam
samtools flagstat $rawbam > $i.stat.raw.txt && \   # 输出比对统计信息
samtools flagstat $sortedbam > $i.stat.q$mq.txt &   

# ******************************************
# 2. 去除 Duplicate Reads
# ******************************************
echo -e "\n\nStep 2: MarkDuplicates"
# 方法1：使用picard 
java -Xmx8g -jar ${EBROOTPICARD}/picard.jar MarkDuplicates  \
    MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=512 VALIDATION_STRINGENCY=LENIENT  \
    INPUT=$sortedbam \
    OUTPUT=$depbam \
    METRICS_FILE=${i}_dedup_metrics.txt && samtools index -@ $nt $depbam
# 方法2：使用GATK 
#gatk  MarkDuplicates -I $sortedbam -O $depbam -M ${i}_dedup_metrics.txt
# 方法3：使用samtools rmdup
#samtools rmdup $sortedbam $depbam

# ******************************************
# 3. Calculate data metrics
# ******************************************
echo -e "\n\nStep 3: Calculate data metrics"
#2.片段inner size，片段选择是否异常
java -Xmx8g -jar ${EBROOTPICARD}/picard.jar CollectInsertSizeMetrics \
      I=$workdir/map/$depbam \
      O=${i}.insert_size_metrics.txt \
      H=${i}.insert_size_histogram.pdf

#3.测序深度与覆盖度统计
singularity exec $sif qualimap bamqc --java-mem-size=8G \
    -nt 2 -nr 10000  \
    -bam $depbam \
    -outfile report.pdf -outformat PDF  \
    -outdir ${i}.depth_and_coverage_report

# ******************************************
# 4. Variant calling
# ******************************************
echo -e "\n\nStep 4: ${i} Variant calling"
gatk  HaplotypeCaller -ERC GVCF -R $REF   \
  -I $depbam \
  --tmp-dir $tmpdir \
  --max-alternate-alleles 4 \
  --native-pair-hmm-threads ${nt} \
  --sample-ploidy 2 \
  -O $outvcf && rm ${rawbam}* ${sortedbam}*

joint call

对于群体gvcf数据，call变异之后还需要合并。可以选择 GATK CombineGVCFs 和 GATK GenotypeGVCFs 或者速度更快的 glnexus

glnexus (GL, genotype likelihood)用于对大规模gvcf 进行 joint call，其由DNAnexus公司开发。

  • 输入：glnexus的输入主要为gvcf文件，支持gatk和deepvariant产生的gvcf文件。
  • 输出：glnexus输出文件为bcf(vcf的二进制格式)，为方便后续使用可以转成vcf文件。

基本使用

# 运行glnexus
module load Singularity/3.7.3
singularity exec $IMAGE/glnexus/1.4.1.sif glnexus_cli -m 500 -t 35 --bed ref.bed --config gatk ../gvcf/*gz > merge_glnexus.bcf

# 将bcf文件转成vcf文件并压缩
module load BCFtools/1.15.1 HTSlib/1.18
bcftools view merge_glnexus.bcf | bgzip -@ 4 -c > all_raw.vcf.gz
  • --dir DIR, -d DIR 过程文件目录，如果不指定，软件默认创建目录为 ./GLnexus.DB，软件运行之前须保证该目录不能存在
  • --config X, -c X 配置文件，默认为 gatk
  • --threads X, -t X 设置线程数，默认为节点所有线程
  • --mem-gbytes X, -m X 预计能使用的内存大小，单位为GB，默认使用节点绝大多数内存容量
  • --bed FILE, -b FILE 3列的bed文件，该区间为软件分析的区间。如果不指定则表示分析所有contig的所有位置
  • --list, -l 指定放置所有gvcf文件路径的文件。如果样本数过多，可以将所有gvcf路径放在一个文件内，每行一个gvcf文件

四 Variant Filtration

过滤流程

1. gatk VariantFiltration

REF=reference.fa
raw_vcf=my_raw.vcf.gz
filter_name=my_snp_filter    # gatk 不会直接过滤，而是根据条件添加上过滤信息。符合则加上PASS标识，否则添加上自定义的filter_name
gatk_vcf=my_gatk.vcf.gz

gatk --java-options "-Xmx16g" VariantFiltration -R $REF -V ${raw_vcf} \
    --filter-expression "QD < 2.0 || FS > 60.0 || MQ < 40.0 "  \
    --cluster-window-size 5 --cluster-size 2 \
    --filter-name ${filter_name} -O ${gatk_vcf}

# 使用 awk 根据gatk结果过滤
zcat ${gatk_vcf} |awk '$0~/#/ || ($7 =="PASS"){print $0}' |gzip - > all.gatked.vcf.gz

注意需要vcf的idx索引文件，如果没有可以自己生成。

gatk IndexFeatureFile -I rwa.vcf.gz

2. vcfutils.pl

vcfutils.pl varFilter -w 5 -W 10 "gzip -dc all.raw.gatked.vcf.gz|" |gzip - >all.varFilter.vcf.gz
# "gzip -dc all.raw.gatked.vcf.gz|" 用于从压缩文件中读取数据

-w 3: indel周围3bp内的SNP会被过滤

-W 10: 10bp内的多个gap会被过滤

3. vcftools

vcftools --gzvcf ${gatked_vcf} --recode --recode-INFO-all \
    --stdout --maf 0.05  --max-missing 0.9 \
    --minDP 4  --maxDP 1000  \
    --not-chr ChrSy --not-chr ChrUn \
    --minQ 30 --minGQ 80 --min-alleles 2  \
    --max-alleles 2 |gzip - > ${vcftools_vcf}
  • --max-missing 0.9：最多允许百分之十的缺失
  • --minQ 30： 质量值
  • --minGQ 80： 最小基因型质量值

Image.png

  • max_missing: 数字越大越严格，0.95表示允许位点在5%的群体中缺失
  • minDP, maxDP： 根据测序深度过滤。
  • min_alleles, max_alleles: 一般不使用多等位位点，需要过滤。
  • clusterWindowSize, clusterSize： 配合使用，如果在五个碱基的窗口内有出现超过两个SNPs，则过滤这些变异。
  • indel_cluster: 基因组10个碱基内有超过两个indel，则过滤。
  • snp_near_indel: 过滤在indel 5个碱基范围内的SNP。

完整流程