首先友情宣傳生信技能樹

• 全國巡講:R基礎(chǔ)，Linux基礎(chǔ)和RNA-seq實戰(zhàn)演練 : 預(yù)告：12月28-30長沙站
• 廣州珠江新城GEO數(shù)據(jù)挖掘滾動開班

題目來源：【生信技能樹】VCF格式文件的shell小練習

首先使用bowtie2軟件自帶的測試數(shù)據(jù)生成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文件區(qū)分成INDEL和SNP條目

$grep '^##' tmp.vcf | grep 'INDEL'
##INFO=<ID=INDEL,Number=0,Type=Flag,Description="Indicates that the variant is an INDEL.">
$grep -v '^#' tmp.vcf | grep 'INDEL' > tmp_INDEL.vcf
$cat tmp_INDEL.vcf | wc -l
69
$grep -v '^#' tmp.vcf | grep -v 'INDEL' > tmp_SNP.vcf
$cat tmp_SNP.vcf | wc -l
36

2. 統(tǒng)計INDEL和SNP條目的各自的平均測序深度

$grep '^##' tmp.vcf | grep 'depth'
##INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth">
$egrep -o 'DP=[0-9]+;' tmp_INDEL.vcf | awk -v FS='=' 'BEGIN{num=0;}{num=num+$2}END{print(num/NR)}'
40.6087
$egrep -o 'DP=[0-9]+;' tmp_SNP.vcf | awk -v FS='=' 'BEGIN{num=0;}{num=num+$2}END{print(num/NR)}'
37

3. 把INDEL條目再區(qū)分成insertion和deletion情況

$awk '$4>$5{print}' tmp_INDEL.vcf > tmp_deletion.vcf
$cat tmp_deletion.vcf |wc -l
68
$awk '$4<$5{print}' tmp_INDEL.vcf > tmp_insertion.vcf
$cat tmp_insertion.vcf | wc -l
1

4. 統(tǒng)計SNP條目的突變組合分布頻率

$awk '{print $4"->"$5}' tmp_SNP.vcf | sort | uniq -c
      7 A->C
      1 A->G
      4 A->T
      2 C->A
      4 C->G
      3 G->A
      2 G->C
      4 G->T
      6 T->A
      1 T->C
      2 T->G

5. 找到基因型不是 1/1 的條目，個數(shù)

$grep -v '^#' tmp.vcf | grep -v '1/1' |wc -l
26
$grep -v '^#' tmp.vcf | grep -v '1/1' |less -SN

6. 篩選測序深度大于20的條目

$egrep -o "DP=[0-9]+;" tmp.vcf |awk '{print(length($0)-4)}'| sort | uniq -c#查看測序深度位數(shù)
      2 1
    103 2
#測序深度最多只有2位數(shù)
$egrep "DP=2[1-9];|DP=[3-9][0-9];" tmp.vcf |wc -l
100
$egrep "DP=2[1-9];|DP=[3-9][0-9];" tmp.vcf |less -SN

7. 篩選變異位點質(zhì)量值大于30的條目

$grep -v '^#' tmp.vcf |awk '$6>30{print}' > tmp_QUALmorethan30.vcf
$cat tmp_QUALmorethan30.vcf | wc -l
99
$less -SN tmp_QUALmorethan30.vcf

8. 組合篩選變異位點質(zhì)量值大于30并且深度大于20的條目

$egrep "DP=2[1-9]|DP=[3-9][0-9]" tmp_QUALmorethan30.vcf | wc -l
97
$egrep "DP=2[1-9]|DP=[3-9][0-9]" tmp_QUALmorethan30.vcf | less -SN

9. 理解DP4=4,7,11,18 這樣的字段，就是 Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases 計算每個變異位點的 AF

1. DP4，高質(zhì)量測序堿基，4個數(shù)字分別對應(yīng)ref-前，ref-后，alt-前，alt-后

• 【理解DP4=4,7,11,18】
  • 參考序列變異位點之前有4個高質(zhì)量測序堿基，之后有7個高質(zhì)量測序堿基;
  • 變異序列變異位點之前有11個高質(zhì)量測序堿基，之后有18個高質(zhì)量測序堿基

2. AF：Allel Frequency 等位基因頻率
   參考：簡書：劉小澤 - VCF格式

AC、AF、AN【和等位基因有關(guān)】： AC：Allele Count該位點變異的等位基因數(shù)目； AF：Allel Frequency 等位基因頻率； AN：Allel Number 等位基因的總數(shù)目
【單看這個不好理解，舉一個二倍體diploid例子：基因型0/1表示為雜合子，該位點只有一個等位基因發(fā)生突變，AF=0.5（在該位點只有50%的等位基因發(fā)生突變），總的等位基因數(shù)目為2；基因型1/1表示為純合子，AC=2，AF=1，AN=2】

所以，
0/0的AF=0；
0/1的AF=0.5；
1/1的AF=1；
另外，AF標簽屬于INFO列

$grep -v "^#" tmp.vcf | awk '{print$10}' |awk -F":" '{print$1}' | sort | uniq -c #查詢樣本基因型
     26 0/1
     79 1/1
$grep -v "^#" tmp.vcf | awk '{if ($0~"1/1"){$8=$8";AF=1";print$0} else if ($0~"0/1"){$8=$8";AF=0.5";print$0}}' > tmp_addAF.vcf

10. 在前面步驟的bam文件里面找到這個vcf文件的某一個突變位點的測序深度表明的那些reads，并且在IGV里面可視化bam和vcf定位到該變異位點。

$grep -v "^#" tmp.vcf |head -10 |less -SN

10.1 如圖，選擇1104突變位點，DP=43，參考序列為C，突變成A

10.2 檢索bam文件
參考：徐洲更hoptop - SAMtools: SAM格式的處理利器
bam文件：
第4列：比對到的位置
第10列：segment序列

$samtools-1.9 view tmp.bam | awk '{if ($3!="*" && $4<=1104 && substr($10,1104-$4+1,1)=="A") print}' | wc -l
45#與vcf文件中DP=43不符合

10.3 tmp.vcf在IGV中可視化

IGV-vcf

10.4 tmp.bam在IGV中可視化

IGV-bam