開(kāi)始做題前需要生成vcf文件先：
# 先生成BAM文件炫隶，用bowtie2
bowtie2 -x ../index/lambda_virus -1 reads_1.fq -2 reads_2.fq | samtools sort -@ 5 -o test2.bam
# 然后生成vcf文件铁追，用bcftools憎兽，當(dāng)然了前提的前提是安裝了bcftools軟件
bcftools mpileup -f ../reference/lambda_virus.fa test2.bam | bcftools call -vm > test2.vcf
想要順利的把題目做完撑蚌，基礎(chǔ)的Linux命令當(dāng)然需要熟練運(yùn)用，更關(guān)鍵的是對(duì) vcf 文件格式和內(nèi)容的了解刽肠。
- 把突變記錄的vcf文件區(qū)分成 INDEL和SNP條目
# 網(wǎng)上有幾種答案溃肪，我比較傾向于這一種。如果是SNP,那么第四列和第五列都應(yīng)該只有一個(gè)字符音五，否則就是InDel
less test2.vcf | grep -v '^#'| awk '{if(length($4)==1 && length($5)==1) print}' | less -S #SNP
less test2.vcf | grep -v '^#'| awk '{if(length($4)!=1 && length($5)!=1) print}' | less -S
# 實(shí)際操作后我發(fā)現(xiàn)上面的代碼也有問(wèn)題惫撰，因?yàn)樯厦鎯尚写a的值分別為36和29下面代碼的行數(shù)是106,理論上兩者之和應(yīng)該等于106，然而明顯不是躺涝。
less test2.vcf | grep -v '^#'| wc
# 經(jīng)過(guò)對(duì)三種情況的對(duì)比后發(fā)現(xiàn)厨钻，第一行代碼沒(méi)問(wèn)題，問(wèn)題在于第二行，邏輯是有問(wèn)題的夯膀，因?yàn)楫?dāng)?shù)?列不等于1時(shí)诗充，第5列可能等于1，且為InDel诱建，所以正確的代碼應(yīng)該是：
less test2.vcf | grep -v '^#'| awk '{if(length($4)!=1) print}' | less -S #INDEL
- 統(tǒng)計(jì)INDEL和SNP條目的各自的平均測(cè)序深度
# 非常樸素的想法：把對(duì)應(yīng)的部分拿出來(lái)然后加和求平均值
less test2.vcf | grep -v '^#'| awk '{if(length($4)!=1) print}' | cut -f 8 | cut -d ";" -f 4 | cut -d "=" -f 2 | awk '{sum+=$1}END{print sum}' #INDEL總測(cè)序深度
less test2.vcf | grep -v '^#'| awk '{if(length($4)==1 && length($5)==1) print}' | cut -f 8 | cut -d ";" -f 1 | cut -d "=" -f 2 | awk '{sum+=$1}END{print sum}' #SNP總測(cè)序深度
- 把INDEL條目再區(qū)分成insertion和deletion情況
# 如果是insert其障，那么第5列比第4列字符數(shù)多；如果是deletion涂佃，那么第5列比第4列字符數(shù)少
less test2.vcf | grep -v '^#'| awk '{if(length($4)!=1) print}' | awk '{if(length($4)<length($5)) print}' | less -S # insert
less test2.vcf | grep -v '^#'| awk '{if(length($4)!=1) print}' | awk '{if(length($4)>length($5)) print}' | less -S # deletion
- 統(tǒng)計(jì)SNP條目的突變組合分布頻率
# 說(shuō)實(shí)話(huà)一開(kāi)始沒(méi)看懂問(wèn)題是什么意思... 意思就是看下SNP的突變類(lèi)型有哪些，各種類(lèi)型的比例是多少蜈敢。
less test2.vcf | grep -v '^#'| awk '{if(length($4)==1 && length($5)==1) print}' | cut -f 4,5 | sort | uniq -c #可以看都有哪些突變類(lèi)型以及每種突變類(lèi)型的數(shù)量
- 找到基因型不是
1/1的條目以及個(gè)數(shù)
# 先找到條目有哪些
less test2.vcf | grep -v "1/1" | grep -v "^#"
# 計(jì)數(shù)
less test2.vcf | grep -v "1/1" | grep -v "^#" | wc
- 篩選測(cè)序深度大于20的條目
# 即DP>20的條目, 反其道而行之辜荠，如下：
less test2.vcf | grep -v "DP=[1-20]" | grep -v "^#" | less -S
- 篩選變異位點(diǎn)質(zhì)量值大于30的條目
# 質(zhì)量值在第6列
less test2.vcf | grep -v "^#" | awk '{if($6>30) print}' | less -S
- 組合篩選變異位點(diǎn)質(zhì)量值大于30并且深度大于20的條目
# 篩選第6列,且滿(mǎn)足第7題條件
less test2.vcf | grep -v "^#" | awk '{if($6>30) print}' | grep -v "DP=[1-20]" | less -S
- 理解
DP4=4,7,11,18這樣的字段，就是 Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases 計(jì)算每個(gè)變異位點(diǎn)的 AF
# 在vcf文件中確實(shí)可以找到DP4這個(gè)tag抓狭，但是到目前為止查詢(xún)的教程都沒(méi)有關(guān)于這個(gè)的解釋?zhuān)瑐X筋伯病。。否过。
- 在前面步驟的bam文件里面找到這個(gè)vcf文件的某一個(gè)突變位點(diǎn)的測(cè)序深度表明的那些reads午笛，并且在IGV里面可視化bam和vcf定位到該變異位點(diǎn)。
其他思考題
- vcf的全稱(chēng)是什么苗桂？是用來(lái)記錄什么信息的標(biāo)準(zhǔn)格式的文本药磺？
? 全稱(chēng)是variant call format，是用來(lái)記錄SNP,INDEL和SV變異信息的文本
- 一般選用哪個(gè)指令查看vcf文件煤伟，為什么不用vim?
? 一般用cat, less就行癌佩；vim可以對(duì)文本進(jìn)行編輯，有誤操作文件的風(fēng)險(xiǎn)
- vcf文件以’##’開(kāi)頭的是什么信息便锨？請(qǐng)認(rèn)真查看這些信息围辙。’#’開(kāi)頭的是什么信息放案？
? 解釋說(shuō)明的信息
- Vcf文件除頭信息姚建，每行有多少列，請(qǐng)?jiān)敿?xì)敘述每行的含義吱殉！請(qǐng)準(zhǔn)確記憶掸冤。
? 每行固定的是8列，后面還有其他列可以添加
? CHROM 和 POS：參考序列名和variant的位置考婴；如果是INDEL的話(huà)贩虾，位置是INDEL的第一個(gè)堿基位置。
? ID：variant的ID沥阱。比如在dbSNP中有該SNP的id缎罢，則會(huì)在此行給出；若沒(méi)有，則用’."表示其為一個(gè)novel variant策精。
? REF 和 ALT：參考序列的堿基 和 Variant的堿基舰始。
? QUAL：Phred格式(Phred_scaled)的質(zhì)量值，表 示在該位點(diǎn)存在variant的可能性咽袜；該值越高丸卷，則variant的可能性越大；計(jì)算方法：
? Phred值 = -10 * log (1-p) p為variant存在的概率; 通過(guò)計(jì)算公式可以看出值為10的表示錯(cuò)誤概率為0.1询刹，該位點(diǎn)為variant的概率為90%谜嫉。
? FILTER：使用上一個(gè)QUAL值來(lái)進(jìn)行過(guò)濾的話(huà)，是不夠的凹联。GATK能使用其它的方法來(lái)進(jìn)行過(guò)濾沐兰，過(guò)濾結(jié)果中通過(guò)則該值為”P(pán)ASS”;
? 若variant不可靠，則該項(xiàng)不為”P(pán)ASS”或”.”蔽挠。
? INFO：這一行是variant的詳細(xì)信息住闯，內(nèi)容很多，以下再具體詳述澳淑。
? FORMAT 和 TTG11B：這兩行合起來(lái)提供了’TTG11B′這個(gè)sample的基因型的信息比原。’TTG11B′代表這該名稱(chēng)的樣品杠巡，是由BAM文件中的@RG下的
? SM 標(biāo)簽決定的量窘。
- 理解format列和樣本列的對(duì)應(yīng)關(guān)系以及GT AD DP的含義。
? GT:樣品的基因型（genotype）,兩個(gè)數(shù)字中間用’/"分 開(kāi)氢拥，這兩個(gè)數(shù)字表示雙倍體的sample的基因型绑改。0 表示樣品中有ref的allele； 1 表示樣品中variant的allele兄一； 2表示有第二個(gè)variant的allele厘线。因此： 0/0 表示sample中該位點(diǎn)為純合的，和ref一致出革； 0/1 表示sample中該位點(diǎn)為雜合的造壮，有ref和variant兩個(gè)基因型； 1/1 表示sample中該位點(diǎn)為純合的骂束，和variant一致耳璧。
? AD:AD(Allele Depth)為sample中每一種allele的reads覆蓋度,在diploid中則是用逗號(hào)分割的兩個(gè)值，前者對(duì)應(yīng)ref基因型展箱，后者對(duì)應(yīng)variant基因型.
? DP:Depth,為sample中該位點(diǎn)的覆蓋度
- Vcf文件第三列如果不是’.’旨枯，出現(xiàn)的rs號(hào)的id是什么？
? dbSNP數(shù)據(jù)庫(kù)的編號(hào)
- Vcf文件的ref混驰，alt列和樣本列的0/1 1/1 或者1/2的聯(lián)系攀隔？
? 可部分參考第5題答案
- Vcf文件一般用什么軟件生成皂贩？請(qǐng)至少說(shuō)出兩個(gè)軟件。請(qǐng)注意不同軟件生成的vcf格式的稍有不同的地方昆汹。
? bcftools軟件明刷、GATK軟件
- Vcf文件一般都比較大，壓縮后的.gz文件用什么指令直接查看而不用解壓后查看满粗？
? 用zcat命令
- 了解gvcf是什么格式辈末，gvcf全稱(chēng)是什么？他與vcf有什么前后聯(lián)系映皆？
? gVCF全稱(chēng)是：Genome Variant Call Format file挤聘。 gVCF和VCF的最大區(qū)別是在于gVCF文件會(huì)記錄所有的點(diǎn)，包括哪些沒(méi)有突變的點(diǎn)捅彻。
- 把a(bǔ)lt列出現(xiàn)’,’的行提取出來(lái)
? 是我沒(méi)理解么檬洞，ALT列不全都是堿基么，怎么會(huì)有‘沟饥，’？
- 請(qǐng)將chrid湾戳、postion贤旷、ref、alt砾脑、format幼驶、樣本列切割出來(lái)生成一個(gè)文本
# 相當(dāng)于提取出第1，2韧衣，4盅藻，5，9畅铭，10列
less test2.vcf | grep -v "^#" | awk '{print $1,$2,$4,$5,$9,$10}' > tmp.txt
- 將一個(gè)含snp,indel信息的vcf拆成一個(gè)只含snp,一個(gè)只含indel信息的2個(gè)vcf文件(可借鑒軟件)
? 參考上面第1題代碼氏淑，直接結(jié)果輸出到一個(gè)新文件即可
- 用指令操作indel的vcf文件，提取indel長(zhǎng)度>4的變異行數(shù)硕噩，存成一個(gè)文本假残。
# 生成indel的vcf文件先
less test2.vcf | grep -v '^#'| awk '{if(length($4)!=1) print}' > InDel.vcf
# 篩選indel長(zhǎng)度大于4的行
cat InDel.vcf | awk '{if(length($5)>4)print}' | less -S
- 用Vcftools過(guò)濾vcf文件，如maf 設(shè)置成0.05炉擅， depth設(shè)置成5-20辉懒，統(tǒng)計(jì)過(guò)濾前后的變異位點(diǎn)的總個(gè)數(shù)
# 如果沒(méi)有需要先下載安裝vcftools軟件，然后--help或者看網(wǎng)上教程學(xué)習(xí)相關(guān)參數(shù)：https://vcftools.github.io/man_latest.html
vcftools --vcf test2.vcf --maf 0.05 --keep-INFO DP=[5-20] --out mtest.vcf
# 運(yùn)行了一下谍失，總是報(bào)錯(cuò)眶俩，我還沒(méi)找到解決辦法暫時(shí)。快鱼。颠印。
- 利用vcftools提取每個(gè)樣本每一個(gè)位點(diǎn)的變異信息和深度信息纲岭，生成一個(gè)矩陣的文件，至少含義以下信息
? 這部分需要學(xué)一點(diǎn)R語(yǔ)言知識(shí)嗽仪，等學(xué)了再來(lái)搞它
