不同样本的WES数据分析时多比对区域是否有差异

Posted on by

做过WES数据分析的小伙伴都知道，最后找到了的很多变异位点，其实都落到了那些富含多比对reads的区域。而这些区域呢，其实是NGS技术的缺点，有些时候可以直接把位于这些区域的突变位点直接硬过滤掉。

对bam文件进行统计，找到那些富含多比对reads的区域

bam文件已经走完了GATK的best practice啦，主要是使用 samtools 和 bedtools 挑选那些被多比对区域的reads，然后根据染色体坐标进行计算覆盖度（测序深度）即可，全部的代码如下：

GENOME=$HOME/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta
samtools faidx $GENOME 
cat $GENOME.fai |grep -v '[-_]' | cut -f1,2 |awk '{print $1"\t1\t"$2}' > hg38.chrom.bed
# 这个 hg38.chrom.bed 文件就是染色体的长度信息啦！
ls N*.bam |while read id;do 
 samtools view $id |\
 perl -alne 'next unless $F[4] eq 0; print if $F[5] =~ /^\d+M$/' |cut -f 3-4|\
 awk '{print $1"\t"$2"\t"$2+150}' | grep -v '[-_]' > ${id%%_*}_bam_qual_0.bed 
 # 为了方便起见，直接把所有的reads长度设置为150bp。
 bedtools coverage -a hg38.chrom.bed -b ${id%%_*}_bam_qual_0.bed -d |\
 perl -lane '{print if $F[4] > 20 }' |cut -f 1,4 > ${id%%_*}_black.list
done

可以看到，每个样品的WES数据分析时多比对区域都是5M左右。通常WES设计的区域是45M，所以这个比例还算是可以接受啦。全部的文件 wc看一下行数：

 5072862 N1_black.list
 4457585 N2_black.list
 4115541 N3_black.list
 4238984 N4_black.list
 4220830 N15_black.list
 4543869 N6_black.list
 4560860 N7_black.list

如果我们把全部的合并起来去冗余后，是 6073693 个位点，其实就是6M区域。也就是说，不同样本的富含多比对reads的区域是类似的， 这个区域具有一定程度的保守性质。

学徒作业：

就选择我4年前的教程：肿瘤全外显子测序数据分析流程大放送，链接是：http://www.bio-info-trainee.com/2735.html

提到的数据集：(SRA) database (accession ID SRA291701)，拿到全部的正常样品的wes数据后走GATK的best practice，使用我前面的代码拿到多比对区域的black.list，探索同样的规律。

历年学徒作业目录如下：

​

Comments are closed.