snp

1. 保留ALT列只突变为一个碱基的行
   1.1 查看snp原文件ALT列的内容

$ grep -v "#" Summit.sorted.markdup.snp.vcf | awk '{print$5}' | sort | uniq -c
 256872 A
    517 *,A
     82 A,*  
     89 A,C
    121 A,G
    112 A,T
 213232 C
    339 *,C
     35 C,*
     45 C,A
     59 C,G
    133 C,T
 213296 G
    290 *,G
     37 G,*
     44 G,A
     28 G,C
    109 G,T
 258061 T
    507 *,T
     89 T,*
     49 T,A
     34 T,C
     32 T,G

1.2 删除ALT列大于一个突变的行
1.21查找不同突变的类型

在nptepad++ 中使用正则表达式查找目标行
.*[AGCT*]{1},[AGCT*]{1}(?!y).*\n
计数2751次匹配
总数等于1.1中非单突变的数量

$ grep -v "#" Summit.sorted.markdup.snp.vcf | awk '{print$5}' | wc -l
944212 #所有突变的行数 不用写NR>1 因为 表头一行前有 #

grep -v "#" Summit.sorted.markdup.snp.vcf | awk  '{if (length($5)==1){print
$0}}'| wc -l
941461  # 突变为单碱基的行数

grep -v "#" Summit.sorted.markdup.snp.vcf | awk  '{if (length($5)==3){print $0}}'| wc -l
2751 #非单突变的行数

1.22 删除非单突变

cat Summit.sorted.markdup.snp.vcf |sed 's/.*[AGCT*]{1},[AGCT*]{1}(?=\t).*\n//' |grep -v "#" | wc -l
944212    #使用sed不能删除非单突变

使用notepad++ 删除非单突变的行后得到过滤后的文件Summit.snp.filtered.vcf
$ grep -v "#" Summit.snp.filtered.vcf | awk '{print$5}' | wc -l
941461 #与原文件单突变数量一样
grep -v "#" Summit.snp.filtered.vcf | awk '{print$5}' | sort | uniq -c
 256872 A
 213232 C
 213296 G
 258061 T #得到的不同碱突变基的类型

1.2.3 #查看GT类型

$ grep -v "#" Summit.snp.filtered.vcf | awk '{print$10}' | awk -F ":" '{print $1}'|sort|uniq -c
 643351 0/1 #杂合突变
 298110 1/1 #纯合突变

2. Summit.snp.filtered.vcf 注释
   2.1 注释前需要构建注释数据库
   下载在对应的注释文件及参考基因组组，
   基因组命名及位置，参考snpeff官方操作说明
   /path/to/snpEff/data/mm37.61/ genes.gtf.gz #注释文件位置，命名为genes.*
   /path/to/snpEff/data/mm37.61/sequences.fa #参考基因组位置，好像只能命名为sequences.fa ，sequences.fna报错
   /path/to/snpEff/data/genomes/mm37.61.fa #参考基因组也可保存在与mm37.61同级的genomes 文件夹下，命名 mm37.61.fa ，mm37.61也是运行snpeff时的参考注释数据名称。
   snpEffectPredictor.bin 为构建数据库生成的文件

cd /path/to/snpEff
java -jar snpEff.jar build -gff3 -v  数据库名称

2.1 注释
Prunusavium 为构建的注释数据库

$ java -Xmx4g -jar /home/Pomgroup/gdp/app/snpeff_v4.4/snpEff/snpEff.jar Prunusavium ./Summit.snp.filtered.vcf > Summit.snp.ann.vcf && echo "well done" || echo "failure"
$ grep -v "#" Summit.snp.filtered.vcf |wc -l
941461
grep -v "#" Summit.snp.filtered.ann.vcf |wc -l
941461  #每行都得到注释

注释成功。
2.2 注释后文件的统计

cat Summit.snp.filtered.ann.genes.txt | awk 'NR > 2{print $4}' | sort | uniq -c
   6050 0
    148 1
     11 2
     10 3
  66985 protein_coding
   1580 pseudogene  
#  在BioType 这一列，为空的不显示内容，内容在两个\t 中间
 cat Summit.snp.filtered.ann.genes.txt | awk -F"\t" 'NR > 2{print $4}' | sort | uniq -c
   6219 （空内容）
  66985 protein_coding
   1580 pseudogene
  指定分隔符 \t

查看snp 注释后的增加的内容

$ sed 's/AC.*;//' Summit.snp.filtered.ann.vcf |grep -v "#" |awk -F "[\t=]" '{print $8}' | sort | uniq -c
 940708 ANN
    348 LOF
    405 NMD
$ sed 's/AC.*;//' Summit.snp.filtered.ann.vcf |grep -v "#" | grep "ANN" |awk -F "[\t=]" '{print $8}' | sort | uniq -c
 940708 ANN

转换颠换情况

$ grep -v "#" Summit.snp.filtered.ann.vcf |awk -F "\t" '{print $4$5}' |sort|uniq -c
  48116 AC
 137168 AG
  56626 AT
  53174 CA
  28661 CG
 148040 CT
 147454 GA
  28239 GC
  53395 GT
  56244 TA
 136877 TC
  47467 TG
             #REF ALT

转换类型为

转换类型
颠换类型