hg19与hg38的差异区域提取

#default

获得了一个新需求：在注释结果中提示突变是否hg19与hg38的差异区域。

需求分析：解读需要更多的辅助信息，提示检出对应突变时，这些突变在hg38中可能是野生型。

方案讨论

Gemini建议的方案是，使用hg19到hg38的liftover对应的chain file，从chain file里获得一致性区域，然后使用hg19的size file来排除这些一致性区域，获得差异区域。

我认为不妥，因为chain file中可能不包含一些小的Gap，也不包含hg19和hg38中的SNP差异。

经过继续讨论，我们确立了差异区域提取方案：

  1. 向探针厂家获得我们使用的WES探针的hg19或hg38版本bed文件；
  2. 使用hg38版本的bed，从hg38的fasta文件中获得对应的exon fasta，提高分析效率。
    • 需要注意，模拟数据使用的wgsim软件，可能无法模拟短于设定长度的序列，因此，需要对bed文件进行调整，确保每条序列大于设定长度。例：我需要模拟PE150数据，就让每个区域大小都大于300bp；
    • 需要对bed进行±50bp处理，以覆盖剪接区域。
  3. 使用wgsim模拟数据，注意设定错误率和突变率都是0；
  4. 将模拟所得的fastq数据比对到hg19，然后进行变异检测，获得SNP和InDel，可以使用hg19的bed±50bp来限制区域；
  5. 所得的SNP和InDel即是基于hg19坐标系的hg19和hg38的在我们的WES探针里的Gap区域；
  6. 对接到VEP注释流程里。

方案实施

处理bed文件

我们的WES探针是安捷伦的，从安捷伦里获得的对应的hg19和hg38的bed文件。

从UCSC获得hg19和hg38的fasta文件。

获得hg38的size file

gunzip hg38.p14.fa
samtools faidx hg38.p14.fa

然后先对hg38的bed进行处理，因为我之后想模拟的是PE150的数据，同时wgsim的默认参数需求插入片段长度（默认值650bp）大于测序长度（300bp）所以这里直接使用±400确保每个区域都大于800bp。这个处理没有关系，因为我们最后是使用hg19的bed来限制输出的SNP和InDel的。

# ±150bp并合并区域
bedtools slop -i S33699751_Covered.bed -b 400 -g hg38.p14.fa.fai > S33699751_Covered.append400.bed

# 排序
bedtools sort -i S33699751_Covered.append400.bed -faidx hg38.p14.fa.fai > S33699751_Covered.append400.sort.bed

# 合并
bedtools merge -i S33699751_Covered.append400.sort.bed > S33699751_Covered.append400.sort.merge.bed

获得模拟数据

接下来从hg38的fasta中，获得我们用于模拟的fasta。

bedtools getfasta -fi hg38.p14.fa -bed S33699751_Covered.append400.sort.merge.bed -fo hg38.select.fa

我需要模拟一个大约30X的数据，需要统计一下大概需要多少reads。

awk '{sum+=$3-$2+1} END {print sum}' S33699751_Covered.append400.sort.merge.bed

如果你用的是和我一样的文件和命令，现在应该会获得同一个值：189077922。因为我想模拟PE150，因此需要将这个值除以300，得到平铺状态下，大概需要630260条reads。但是显然，这个数值是不对的，因为我们有多个区域，因此我就直接将预估的reads数改为700000，在30X下，即21000000。这个21000000是计算单端的，为了容错，狠心一点，翻个倍，改为42000000😄，这样预估出来的平均深度应该是60X的。

其他参数，咱也不懂就不要改。

wgsim -e 0 -r 0 -R 0 -d 500 -s 50 -N 42000000 -1 150 -2 150 hg38.select.fa out_1.fq out_2.fq

# 压缩结果文件
gzip out_1.fq 
gzip out_2.fq

变异检测

这一步是比对到hg19和进行变异检测，就不细说了，大家都有自己现成的方案。我们的变异检测流程里，是已经对exon区域进行了±50bp处理的，因此结果vcf文件其实可以直接使用。

处理VCF

VEP可以直接把VCF文件作为注释源。我需要给这个VCF文件的INFO列增加标记信息，并且删除FORMAT列和样本列，节省空间。

首先我们移除FORMAT列和样本信息列

bcftools view -Ov -s - input.vcf --force-samples > no_sample.vcf

接下来增加INFO列的标记

awk '
BEGIN { OFS="\t" } # 设置输出字段分隔符为 Tab

# 处理VCF头部行
/^##/ {
    print $0 # 打印所有以 "##" 开头的行
    next     # 跳到下一行
}

# 处理VCF列名行（#CHROM ...）
/^#CHROM/ {
    # 在 #CHROM 行之前插入新的 ##INFO 定义
    print "##INFO=<ID=HG38,Number=1,Type=String,Description=\"hg38 wild type mark.\">"
    print $0 # 打印原始的 #CHROM 行
    next     # 跳到下一行
}

# 处理数据行
{
    # 给INFO列增加信息
    $8 = $8 ";HG38=WT"

    # 重新打印前8列
    print $1,$2,$3,$4,$5,$6,$7,$8
}
' no_sample.vcf > hg19_hg38wt.vcf

使用bgzip压缩VCF并建立索引

bgzip hg19_hg38wt.vcf
tabix -p vcf hg19_hg38wt.vcf.gz

VEP注释

可以使用VEP将这个数据库注释到结果中，下面是示例

vep \
	-i input.vcf \
	-o output.vcf \
	----custom file=hg19_hg38wt.vcf.gz,short_name=Diff,format=vcf,type=exact,coords=0,fields=HG38 \
	--fields Diff_HG38

这样，在我们基于hg19的WES分析流程里，就可以在结果中看到一个突变在hg38中是野生型的提示。

请注意

这个新增注释是hg19的错误注释的充分不必要条件。hg19的转录本中可能会把一些hg38的内含子区划分到外显子中，因此会造成转录本水平的差异，但在基因组水平可能不存在差异。