VCF格式文件的shell小练习

首先使用bowtie2软件自带的测试数据生成sam/bam文件，还有vcf文件
代码如下：

mkdir -p ~/biosoft
cd ~/biosoft
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh 
bash Miniconda3-latest-Linux-x86_64.sh 
source ~/.bashrc 
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes
conda create -y -n test
conda activate test
conda install -y samtools bcftools

wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.4.3/bowtie2-2.3.4.3-linux-x86_64.zip 
unzip bowtie2-2.3.4.3-linux-x86_64.zip 
cd ~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads

../../bowtie2 -x ../index/lambda_virus -1 reads_1.fq -2 reads_2.fq | samtools sort -@ 5 -o tmp.bam
bcftools mpileup -f ../reference/lambda_virus.fa tmp.bam | bcftools call -vm > tmp.vcf
LINUX练习题
  1. 把突变记录的vcf文件区分成 INDEL和SNP条目
grep 'INDEL' tmp.vcf | grep -v '#'   # INDEL条目
grep -v -e 'INDEL' -v -e '#' tmp.vcf   # SNP条目
  1. 统计INDEL和SNP条目的各自的平均测序深度
cat tmp.vcf | grep -v "^#" | grep 'INDEL' | cut -f 8 | cut -d ";" -f 4 | awk '{print substr($0, 4)}' | awk '{sum += $0} END {print sum/NR}'   # INDEL条目
cat tmp.vcf | grep -v "^#" | grep -v 'INDEL' | cut -f 8 | cut -d ";" -f 1 | awk '{print substr($0, 4)}' | awk '{sum += $0} END {print sum/NR}'   # SNP条目
  1. 把INDEL条目再区分成insertion和deletion情况
grep -v '^#' tmp.vcf | grep 'INDEL' | awk '{if ($4 > $5) print $0}' | less -SN   # Deletion
grep -v '^#' tmp.vcf | grep 'INDEL' | awk '{if ($4 < $5) print $0}' | less -SN   # Insertion
  1. 统计SNP条目的突变组合分布频率
grep -v '^#' tmp.vcf | grep -v 'INDEL' | cut -f 4-5 | tr '\t' '-' | sort | uniq -c
# grep -v '^#' tmp.vcf | grep -v 'INDEL' | awk '{print $4"-"$5}'| sort | uniq -c
  1. 找到基因型不是 1/1 的条目，个数
grep -v '^#' tmp.vcf | awk '{if ($10 !~ "1/1") print}'
grep -v '^#' tmp.vcf | awk '{if ($10 !~ "1/1") print}' | wc -l
  1. 筛选测序深度大于20的条目
egrep 'DP=2[1-9]|DP=[3-9][0-9]' tmp.vcf | less -SN
  1. 筛选变异位点质量值大于30的条目
grep -v '^#' tmp.vcf | awk '{if ($6 > 30) print}' | less -SN
  1. 组合筛选变异位点质量值大于30并且深度大于20的条目
grep -v '^#' tmp.vcf | awk '{if ($6 > 30) print}' | egrep 'DP=2[1-9]|DP=[3-9][0-9]' | less -SN
  1. 理解DP4=4,7,11,18 这样的字段，就是 Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases 计算每个变异位点的 AF
grep -v "^#" tmp.vcf | cut -f 8 | egrep -o 'DP4=[0-9]+,[0-9]+,[0-9]+,[0-9]+' | egrep -o '[0-9]+,[0-9]+,[0-9]+,[0-9]+' | awk -v FS="," '{print($3+$4)/($1+$2+$3+$4)}'
  1. 在前面步骤的bam文件里面找到这个vcf文件的某一个突变位点的测序深度表明的那些reads，并且在IGV里面可视化bam和vcf定位到该变异位点。
less -SN tmp.vcf

截图前几行：

以34行和35行为例：

在BAM文件中定位reads:

samtools view tmp.bam | awk '{if ($6 != "*" && $4 <= 1104 && substr($10, 1104-$4+1, 1) == "A") print}' | less -SN
samtools view tmp.bam | awk '{if ($6 != "*" && $4 <= 1344 && substr($10, 1344-$4+1, 1) == "T") print}' | less -SN

IGV可视化：

其它思考题
  1. vcf的全称是什么？是用来记录什么信息的标准格式的文本？
    Variant Call Format。记录基因序列变异信息的文本文件。
  2. 一般选用哪个指令查看vcf文件，为什么不用vim?
    less -SN，不用vim可能是处理大文件时效率低。
  3. vcf文件以’##’开头的是什么信息？请认真查看这些信息。’#’开头的是什么信息？
    以 “##” 开头的是注释信息（meta-information）。以 “#” 开头的是标题行。
  4. Vcf文件除头信息，每行有多少列，请详细叙述每行的含义！请准确记忆。
    每行有8个固定的、必需的列。第9列开始是基因型信息。
  1. 理解format列和样本列的对应关系以及GT AD DP的含义。
    format列和样本列的数据是对应的，前者为格式，后者为格式对应的数据。冒号分隔。
    GT: 样本基因型。两个数字中间用 “/” 分开，这两个数字表示二倍体sample的基因型。0 表示样品中有ref的allele；1 表示样品中variant的allele；2表示有第二个variant的allele。0/0 表示sample中该位点为纯合的，和ref一致；0/1 表示sample中该位点为杂合的，有ref和variant两种基因型；1/1 表示sample中该位点为纯合的，和variant一致。
    AD: Allele Depth，为sample中每一种allele的reads覆盖度，在二倍体中是用逗号分割的两个值，前者对应ref基因型，后者对应variant基因型。
    DP: sample中该位点的覆盖度。
  2. Vcf文件第三列如果不是’.’，出现的rs号的id是什么？
    dbSNP数据库中该SNP的ID号。
  3. Vcf文件的ref，alt列和样本列的0/1 1/1 或者1/2的联系？
  4. Vcf文件一般用什么软件生成？请至少说出两个软件。请注意不同软件生成的vcf格式的稍有不同的地方。
    bcftools，GATK4。
  5. Vcf文件一般都比较大，压缩后的.gz文件用什么指令直接查看而不用解压后查看？
    zcat命令 
gzip tmp.vcf
zcat tmp.vcf.gz
gunzip tmp.vcf.gz
  1. 了解gvcf是什么格式，gvcf全称是什么？他与vcf有什么前后联系？
    Genome Variant Call Format。GVCF文件包含所有位点的记录，既包括存在变异的位点，也包括不存在变异的位点。这样做的目的是为了方便后续的群体分析。
  2. 把alt列出现’,’的行提取出来
  3. 请将chrid、postion、ref、alt、format、样本列切割出来生成一个文本
grep -v "^##" tmp.vcf | cut -f 1,2,4,5,9- > tmp_cut.txt
  1. 将一个含snp,indel信息的vcf拆成一个只含snp,一个只含indel信息的2个vcf文件。可借鉴软件
cat tmp.vcf | sed '/INDEL/d' > tmp_SNP.vcf   # SNP
cat tmp.vcf | awk '/#|INDEL/{print $0}' > tmp_INDEL.vcf   # INDEL

bcftools view tmp.vcf -O u -o tmp.bcf
bcftools view tmp.bcf
bcftools view -v snps tmp.bcf | less -SN   # SNP
bcftools view -v indels tmp.bcf | less -SN   # INDEL

vcftools --remove-indels --recode --recode-INFO-all --vcf tmp.vcf --stdout > tmp.snp.vcf   # SNP
vcftools --keep-only-indels --recode --recode-INFO-all --vcf tmp.vcf --stdout > tmp.indel.vcf   # INDEL
  1. 用指令操作indel的vcf文件，提取indel长度>4的变异行数，存成一个文本。
grep -v "^#" tmp.vcf | awk '{if (length($5) > 4) print}'  > 4.txt
  1. 用Vcftools过滤vcf文件，如maf 设置成0.05， depth设置成5-20，统计过滤前后的变异位点的总个数
vcftools --recode \
--recode-INFO-all \
--vcf tmp.vcf \
--maf 0.05 \
--minDP 5 --maxDP 20 \
--stdout > filter.vcf
# 会报错，需要再回顾学习
  1. 利用vcftools提取每个样本每一个位点的变异信息和深度信息，生成一个矩阵的文件，至少含义以下信息

Eg:

CHRID POSTION SAMPLE_DP SAMPLE_GT
chr1 1010 5 AA
  1. 提取出变异位点上样本有纯和突变的行数
  2. 统计一下1号染色体上的变异总个数。
  3. 提取一下BRCA1基因上发生变异的行数，如果是人类wes变异结果文件
  4. 统计vcf文件各样本的缺失率，如果是多个样本的群体call结果。
    进阶思考题：可视化vcf中变异位点质量值，横坐标是质量值，纵坐标是百分比或者该质量值时的变异位点的总个数，可以化成线图或者点图（可能需要学一些代码，还有R画图）

参考：
http://www.bio-info-trainee.com/3577.html
https://samtools.github.io/hts-specs/VCFv4.2.pdf
https://www.jianshu.com/p/eadfcd7e23ee
http://seqanswers.com/forums/showthread.php?t=10181
http://samtools.sourceforge.net/mpileup.shtml
https://cloud.tencent.com/developer/article/1054971

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
平台声明：文章内容（如有图片或视频亦包括在内）由作者上传并发布，文章内容仅代表作者本人观点，简书系信息发布平台，仅提供信息存储服务。

推荐阅读更多精彩内容