pchi縮寫(xiě)是什么意思

PCHI英文含義

1、PCHI的英文全稱(chēng)：Permanent Childhood Hearing Impairment | 中文意思：───兒童持續性聽(tīng)損傷

2、PCHI的英文全稱(chēng)：Permanent Congenital Hearing Impairment | 中文意思：───永久性先天性聽(tīng)力障礙

3、PCHI的英文全稱(chēng)：Partners Community Healthcare, Inc. (est. 1994) | 中文意思：───合作伙伴社區醫療保健公司 (est。 1994)

4、PCHI的英文全稱(chēng)：Philippine Chamber of Handicraft Industries | 中文意思：───菲律賓手工業(yè)的分庭

VCF文件參數解讀

    VCF文件的開(kāi)頭是整體注釋信息，通常以##作為起始，其后一般接以FILTER，INFO，FORMAT等字樣。

    例如：以##FILTER開(kāi)頭的行，表示注釋VCF文件當中第7列中縮寫(xiě)詞的說(shuō)明；##INFO開(kāi)頭的行注釋VCF第8列中的縮寫(xiě)字母說(shuō)明，比如AF代表Allele Frequency也就是等位基因頻率；##FORMAT開(kāi)頭的行注釋VCF第9列中的縮寫(xiě)字母說(shuō)明；另外還有其他的一些信息，文件版本"fileformat=VCFv4.0"等等。還能看到一些歷史命令，通過(guò)這些命令可以知道這個(gè)vcf文件是如何得到的。

    各列之間用tab空白隔開(kāi)；前面9列為固定列，第10列開(kāi)始為樣品信息列，可以無(wú)限多個(gè)；圖示樣品信息列有130個(gè)

    CHROM    POS    ID    REF    ALT    QUAL    FILTER    INFO    FORMAT

    后面的列都為樣品基因型信息列

     1.CHROM     記錄染色體編號

     2.POS    記錄變異位點(diǎn)在參考基因組中的位置。如果是SNP的話(huà)，POS即SNP的位置；如果是INDEL的話(huà)，位置是INDEL的第一個(gè)堿基位置。

     3.ID     SNP/INDEL的ID, 如在dbSNP中有該SNP的id，則會(huì )在此行給出；若沒(méi)有，則用’.'表示其為一個(gè)novel variant 新變異,dbSNP編號通常以rs開(kāi)頭，一般只有人類(lèi)基因組才有dbSNP編號

     INDEL 指的是在基因組的某個(gè)位置上所發(fā)生的small deletion,small inverion小片段序列的**入或者刪除，其長(cháng)度通常在50bp以下

     4.REF     參考基因組該位置堿基類(lèi)型，必須是A,C,G,T,N    N表示不確定堿基，SNP應該一個(gè)位點(diǎn)就是一個(gè)堿基

     5.ALT     與參考序列比較，發(fā)生突變的變異堿基類(lèi)型，必須是A,C,G,T,N，.    多個(gè)用逗號分割。"." 表示這個(gè)地方?jīng)]有reads覆蓋為缺失。

     6.QUAL      變異位點(diǎn)檢測質(zhì)量值，越高越可靠。表示在該位點(diǎn)存在variant的可能性，該值越高，則variant的可能性越大 

        等于-10*log10(該變異位點(diǎn)檢測錯誤的概率)。 用 . 表示，是質(zhì)量值沒(méi)有輸出，不代表質(zhì)量值為0

           log0.1表示10的多少次方等于0.1，即為-1；10的-1次方為十分之一，10的-2次方為一百分之一

     7.FILTER     如果該位點(diǎn)通過(guò)過(guò)濾標準那么我們可以在該列標記為"PASS",說(shuō)明該列質(zhì)量值高。

     8.  INFO為variant的詳細信息     字段的意思可以在header里搜索去看

     上面vcf 中INFO全為“.”了，是因為用 vcftools 某步過(guò)濾SNP輸出文件時(shí)用了 --recode ，這樣就不輸出info信息，以 . 代替了，想輸出info，可以--recode-INFO xx（如MQ）  或者 --recode-INFO-all （所有info全部輸出）

    #DP-read depth：樣本在這個(gè)位置的reads覆蓋度。是一些reads被過(guò)濾掉后的覆蓋度。DP4:高質(zhì)量測序堿基，位于REF或者ALT前后

    #QD：通過(guò)深度來(lái)評估一個(gè)變異的可信度。Variant call confidence normalized by depth of sample reads supporting a variant       

    #MQ：表示覆蓋序列質(zhì)量的均方值RMS Mapping Quality

    #FQ：phred值關(guān)于所有樣本相似的可能性

     #AC，AF 和 AN：AC(Allele Count) 表示該Allele的數目；AF(Allele Frequency) 表示Allele的頻率； AN(Allele Number) 表示Allele的總數目。

     對于1個(gè)diploid sample（雙倍體）而言：則基因型 0/1 表示sample為雜合子，Allele數為1 (雙倍體的sample在該位點(diǎn)只有1個(gè)等位基因發(fā)生了突變)，Allele的頻率為0.5 (雙倍體的sample在該位點(diǎn)只有50%的等位基因發(fā)生了突變)，總的Allele為2；                                                                    基因型 1/1 則表示sample為純合的，Allele數為2，Allele的頻率為1，總的Allele為2。

    #MLEAC：Maximum likelihood expectation (MLE) for the allele counts (not necessarily the same as the AC), for each ALT allele, in the same order as listed

    #MLEAF：Maximum likelihood expectation (MLE) for the allele frequency (not necessarily the same as the AF), for each ALT allele, in the same order as listed

    #BaseQRankSum  比較支持變異的堿基和支持參考基因組的堿基的質(zhì)量，負值表示支持變異的堿基質(zhì)量值不及支持參考基因組的，

     正值則相反，支持變異的質(zhì)量值好于參考基因組的。0表示兩者無(wú)明顯差異。

    #FS  使用F檢驗來(lái)檢驗測序是否存在鏈偏好性。鏈偏好性可能會(huì )導致變異等位基因檢測出現錯誤。輸出值Phred-scaled p-value，值越大越可能出現鏈偏好性。

    #InbreedingCoeff    使用似然法檢驗樣本間的近交系數（又或者稱(chēng)為近親關(guān)系）。值越高越可能是近親繁殖。

    #MQRankSum  比較支持變異的序列和支持參考基因組的序列的質(zhì)量，負值表示支持變異的堿基質(zhì)量值不及支持參考基因組的，只針對雜合。

     正值則相反，支持變異的質(zhì)量值好于參考基因組的。0表示兩者無(wú)明顯差異。實(shí)際應用中一般過(guò)濾掉較小的負值。

    #BaseCounts  所有樣本在變異位點(diǎn)ATCG的數量

    #ClippingRankSum  同前面兩個(gè)類(lèi)似，負值表示支持變異的read有更的的hard-clip堿基，正值表示支持參考基因組的的read有更多的hard-clip。0最好，無(wú)論是正值還是負值都表示可能可能存在人為偏差。

    #ReadPosRankSum    檢測變異位點(diǎn)是否有位置偏好性（是否存在于序列末端，此時(shí)往往容易出錯）。最佳值為0，表示變異與其在序列上的位置無(wú)關(guān)。負值表示變異位點(diǎn)更容易在末端出現，正值表示參考基因組中的等位基因更容易在末端出現。

    #ExcessHet  檢測這些樣本的相關(guān)性，與InbreedingCoeff相似，值越大越可能是錯誤。

    #LikelihoodRankSum  評價(jià)支持變異和ref的序列與best hyplotype的匹配性，0為最佳值。負值表示支持變異的read匹配度不及支持ref的匹配度，正值則相反。值越大表示越可能是出現了錯誤。

    #HaplotypeScore    分數越高越可能出現錯誤。Higher scores are indicative of regions with bad alignments, typically leading to artifactual SNP and indel calls.

    #SOR：也是一個(gè)用來(lái)評估是否存在鏈偏向性的參數，相當于FS的升級版。The StrandOddsRatio annotation is one of several methods that aims to evaluate whether there is strand bias in the data. It is an updated form of the Fisher Strand Test that is better at taking into account large amounts of data in high coverage situations. It is used to determine if there is strand bias between forward and reverse strands for the reference or alternate allele. The reported value is ln-scaled.

    #IS：**入缺失或部分**入缺失的reads允許的最大數量

    #G3：ML 評估基因型出現的頻率

    #HWE：chi^2基于HWE的測試p值和G3

    #CLR：在受到或者不受限制的情況下基因型出現可能性log值

    #UGT：最可能不受限制的三種基因型結構

    #CGT：最可能受限制三種基因型的結構

    #PV4：四種P值的誤差，分別是（strand、baseQ、mapQ、tail distance bias）

    #INDEL：表示該位置的變異是**入缺失

    #PC2：非參考等位基因的phred（變異的可能性）值在兩個(gè)分組中大小不同

    #PCHI2：后加權chi^2，根據p值來(lái)測試兩組樣本之間的聯(lián)系

    #QCHI2：Phred scaled PCHI2

    #PR：置換產(chǎn)生的一個(gè)較小的PCHI2

    #QBD：Quality by Depth，測序深度對質(zhì)量的影響

    #RPB：序列的誤差位置（Read Position Bias）

    #MDV：樣本中高質(zhì)量非參考序列的最大數目

    #VDB：Variant Distance Bias，RNA序列中過(guò)濾人工拼接序列的變異誤差范圍

     9.FORMAT 為后面10列信息的說(shuō)明列，通常以" ："隔開(kāi)各個(gè)縮寫(xiě)詞。

     10 列（包含）以后 為樣品基因型列，各信息以"："分隔與FORMAT列一一對應；

（不確定 1/0與0/1 , 1/2與2/1 , 2/3與3/2 是否為一個(gè)意思，猜測可能是一個(gè)意思，沒(méi)有去深究）

    在過(guò)濾后只剩SNP的vcf文件中，GT只會(huì )存在 0/0    0/1      1/1       0（參考基因組等位基因類(lèi)型）和1（樣品的一種變異等位基因類(lèi)型）

     像下圖，還存在除SNP外其他類(lèi)型的變異，所以GT存在1/2，2/2等

     AD 和DP： AD(Allele Depth)為sample中在此位置支持每種堿基型的reads深度，用逗號分割，前者對應ref基因型，后者對應variant基因型； DP（Depth）為sample中該位點(diǎn)的覆蓋度，為該變異位點(diǎn)的深度和，也就是AD兩個(gè)數字的和。

     GQ ： 基因型質(zhì)量值        Phred值 = -10 * log (p)  p為基因型錯誤的概率  越高越可靠

     PL ： 指定的三種基因型的似然值。這三種指定的基因型為(0/0,0/1,1/1)，這三種基因型的概率總和為1。數值越小代表基因型越可靠，最小的數字對應的基因型判讀為該樣品的最可能的基因型。比如最后一列285,0,105，分別對應基因型0/0,0/1,1/1，說(shuō)明0/1為可能的基因型。

     PGT   PID  也看了，沒(méi)咋懂，不記錄了

參考：

https://www.jianshu.com/p/1726696e54e5

https://www.jianshu.com/p/13f162636164

https://www.omicsclass.com/article/847

https://www.jianshu.com/p/bf0d27368eb9

https://www.jianshu.com/p/13f162636164

samtools使用大全

samtools是一個(gè)用于操作sam和bam文件（通常是短序列比對工具如bwa,bowtie2,hisat2,tophat2等等產(chǎn)生的，具體格式可以在消息框輸入“SAM”查看）的工具合集，包含有許多命令。以下是常用命令的介紹。

1.View

view命令的主要功能是：將sam文件與bam文件互換；然后對bam文件進(jìn)行各種操作，比如數據的排序(sort)和提取(這些操作 是對bam文件進(jìn)行的，因而當輸入為sam文件的時(shí)候，不能進(jìn)行該操作)；最后將排序或提取得到的數據輸出為bam或sam（默認的）格式。

bam文件優(yōu)點(diǎn)：bam文件為二進(jìn)制文件，占用的磁盤(pán)空間比sam文本文件??；利用bam二進(jìn)制文件的運算速度快。

view命令中，對sam文件頭部（序列ID）的輸入(-t或-T）和輸出(-h)是單獨的一些參數來(lái)控制的。

Usage: samtools view [options] | [region1 [...]]

默認情況下不加 region，則是輸出所有的 region.

options:

-b output BAM

  默認下輸出是 SAM 格式文件，該參數設置輸出 BAM 格式

-h print header for the SAM output

  默認下輸出的 sam 格式文件不帶 header，該參數設定輸出sam文件時(shí)帶 header 信息

-H print header only (no alignments)

  僅僅輸出文件的頭

-S input is SAM

  默認下輸入是 BAM 文件，若是輸入是 SAM 文件，則最好加該參數，否則有時(shí)候會(huì )報錯。

-u uncompressed BAM output (force -b)

  該參數的使用需要有-b參數，能節約時(shí)間，但是需要更多磁盤(pán)空間。

-c Instead of printing the alignments, only count them and print the

  total number. All filter options, such as ‘-f’, ‘-F’ and ‘-q’ ,  are taken into account.

  過(guò)濾功統計功能

-c print only the count of matching records

-L FILE  output alignments overlapping the input BED FILE [null]

-t FILE  list of reference names and lengths (force -S) [null]

  使用一個(gè)list文件來(lái)作為header的輸入

-T FILE  reference sequence file (force -S) [null]

  使用序列fasta文件作為header的輸入

-o FILE  output file name [stdout]

-F INT  filtering flag, 0 for unset [0]

  Skip alignments with bits present in INT [0]

  數字4代表該序列沒(méi)有比對到參考序列上

  數字8代表該序列的mate序列沒(méi)有比對到參考序列上

  過(guò)濾功能。如F12過(guò)濾只有雙端map的

-q INT  minimum mapping quality [0]

    比對的最低質(zhì)量值，一般認為20就為unique比對了，可以結合上述-bF參數使用使用提取特定的比對結果

例子：

將sam文件轉換成bam文件

samtools view -bS abc.sam >abc.bam

BAM轉換為SAM

samtools view -h -o out.sam out.bam

提取比對到參考序列上的比對結果

samtools view -bF 4 abc.bam >abc.F.bam

提取paired reads中兩條reads都比對到參考序列上的比對結果，只需要把兩個(gè)4+8的值12作為過(guò)濾參數即可

samtools view -bF 12 abc.bam >abc.F12.bam

提取沒(méi)有比對到參考序列上的比對結果

samtools view -bf 4 abc.bam >abc.f.bam

提取bam文件中比對到caffold1上的比對結果，并保存到sam文件格式

samtools view abc.bam scaffold1 >scaffold1.sam

提取scaffold1上能比對到30k到100k區域的比對結果

samtools view abc.bam scaffold1:30000-100000 $gt; scaffold1_30k-100k.sam

根據fasta文件，將 header 加入到 sam 或 bam 文件中

samtools view -T genome.fasta -h scaffold1.sam >scaffold1.h.sam

2.Sort

sort對bam文件進(jìn)行排序。一些軟件需要sort的bam或者sam文件，如stringtie，所以必須要sort使用；求depth時(shí)，也必須要sort；

Usage: samtools sort [-n] [-m ] 

-m 內存參數默認下是 500,000,000 即500M（不支持K，M，G等縮寫(xiě)）。對于處理大數據時(shí)，如果內存夠用，則設置大點(diǎn)的值，以節約時(shí)間。

-n 設定排序方式按short reads的ID排序。默認下是按序列在fasta文件中的順序（即header）和序列從左往右的位點(diǎn)排序。

例子：

samtools sort accept.bam accept.sort最終產(chǎn)生accept.sort.bam

3.merge

將2個(gè)或2個(gè)以上的已經(jīng)sort了的bam文件融合成一個(gè)bam文件。融合后的文件已經(jīng)sort過(guò)了的。

Usage:  samtools merge [-nr] [-h inh.sam] [...]

Options: -n      sort by read names

        -r      attach RG tag (inferred from file names)

        -u      uncompressed BAM output

        -f      overwrite the output BAM if exist

        -1      compress level 1

        -R STR  merge file in the specified region STR [all]

        -h FILE  copy the header in FILE to [in1.bam]

例子：

4.index

必須對bam文件進(jìn)行默認情況下的排序后，才能進(jìn)行index。否則會(huì )報錯。

建立索引后將產(chǎn)生后綴為.bai的文件，用于快速的隨機處理。很多情況下需要有bai文件的存在，特別是顯示序列比對情況下。比如samtool的tview命令就需要；gbrowse2顯示reads的比對圖形的時(shí)候也需要。IGV顯示比對情況也需要。

Usage: samtools index [out.index]

例子：

以下兩種命令結果一樣

$ samtools index abc.sort.bam

$ samtools index abc.sort.bam abc.sort.bam.bai

5.faidx

對fasta文件建立索引,生成的索引文件以.fai后綴結尾。該命令也能依據索引文件快速提取fasta文件中的某一條（子）序列

Usage: samtools faidx [ [...]]

對基因組文件建立索引，方便提取序列。

例子：$ samtools faidx genome.fasta

由于有索引文件，可以使用以下命令很快從基因組中提取到fasta格式的子序列

$ samtools faidx genome.fasta scffold_10 >scaffold_10.fasta

6.tview

tview能直觀(guān)的顯示出reads比對基因組的情況，和基因組瀏覽器有點(diǎn)類(lèi)似。

需要事先利用利用上面講的sort和建index命令執行完后，用下述命令。

Usage: samtools tview [ref.fasta]

出參考基因組的時(shí)候，會(huì )在第一排顯示參考基因組的序列，否則，第一排全用N表示。

按下 g ，則提示輸入要到達基因組的某一個(gè)位點(diǎn)。例子“scaffold_10:1000"表示到達第

10號scaffold的第1000個(gè)堿基位點(diǎn)處。

使用H(左）J（上）K（下）L（右）移動(dòng)顯示界面。大寫(xiě)字母移動(dòng)快，小寫(xiě)字母移動(dòng)慢。

使用空格建向左快速移動(dòng)（和 L 類(lèi)似），使用Backspace鍵向左快速移動(dòng)（和 H 類(lèi)似）。

Ctrl+H 向左移動(dòng)1kb堿基距離； Ctrl+L 向右移動(dòng)1kb堿基距離

可以用顏色標注比對質(zhì)量，堿基質(zhì)量，核苷酸等。30～40的堿基質(zhì)量或比對質(zhì)量使用白色表示；

20～30**；10～20綠色；0～10藍色。

使用點(diǎn)號'.'切換顯示堿基和點(diǎn)號；使用r切換顯示read name等

還有很多其它的使用說(shuō)明，具體按 ？ 鍵來(lái)查看。

7.flagstat

給出BAM文件的比對結果

Usage: samtools flagstat

$ samtools flagstat example.bam

11945742 + 0 in total (QC-passed reads + QC-failed reads)

#總共的reads數

0 + 0 duplicates

7536364 + 0 mapped (63.09%:-nan%)

#總體上reads的匹配率

11945742 + 0 paired in sequencing

#有多少reads是屬于paired reads

5972871 + 0 read1

#reads1中的reads數

5972871 + 0 read2

#reads2中的reads數

6412042 + 0 properly paired (53.68%:-nan%)

#完美匹配的reads數：比對到同一條參考序列，并且兩條reads之間的距離符合設置的閾值

6899708 + 0 with itself and mate mapped

#paired reads中兩條都比對到參考序列上的reads數

636656 + 0 singletons (5.33%:-nan%)

#單獨一條匹配到參考序列上的reads數，和上一個(gè)相加，則是總的匹配上的reads數。

469868 + 0 with mate mapped to a different chr

#paired reads中兩條分別比對到兩條不同的參考序列的reads數

243047 + 0 with mate mapped to a different chr (mapQ>=5)

#同上一個(gè)，只是其中比對質(zhì)量>=5的reads的數量

8.depth

得到每個(gè)堿基位點(diǎn)的測序深度,并輸出到標準輸出，所以要用大于號追加到一個(gè)文件。

Usage: bam2depth [-r reg] [-q baseQthres] [-Q mapQthres] [-b in.bed] [...]

-r 后面跟染色體號（region）

-q ：計算深度時(shí)要求測序堿基質(zhì)量最低質(zhì)量值

-Q ：計算深度時(shí)要求比對的最低質(zhì)量值

注意：做depth之前必須做samtools index；

例子

samtools depth accept.bam >depth

9.其他命令

reheader：替換bam文件的頭

$ samtools reheader

idxstats ：統計一個(gè)表格，4列，分別為”序列名，序列長(cháng)度，比對上的reads數，unmapped reads number。第4列應該是paired reads中有一端能匹配到該scaffold上，而另外一端不匹配到任何scaffolds上的reads數。

$ samtools idxstats

rmdup：將由PCR duplicates獲得的reads去掉，并只保留最高比對質(zhì)量的read。

Usage:  samtools rmdup [-sS] 

-s 對single-end reads。默認情況下，只對paired-end reads

-S 將Paired-end reads作為single-end reads處理。

10. 將bam文件轉換為fastq文件

有時(shí)候，我們需要提取出比對到一段參考序列的reads，進(jìn)行小范圍的分析，以利于debug等。這時(shí)需要將bam或sam文件轉換為fastq格式。

該網(wǎng)站提供了一個(gè)bam轉換為fastq的程序：http://www.hudsonalpha.org/gsl/information/software/bam2fastq

$ wget http://www.hudsonalpha.org/gsl/static/software/bam2fastq-1.1.0.tgz

$ tar zxf bam2fastq-1.1.0.tgz

$ cd bam2fastq-1.1.0

$ make

$ ./bam2fastq

11. mpileup

samtools還有個(gè)非常重要的命令mpileup，以前為pileup。該命令用于生成bcf文件，再使用bcftools進(jìn)行SNP和Indel的分析。bcftools是samtool中附帶的軟件，在samtools的安裝文件夾中可以找到。

最常用的參數有2：

-f 來(lái)輸入有索引文件的fasta參考序列； -g 輸出到bcf格式。用法和最簡(jiǎn)單的例子如下

Usage: samtools mpileup [-EBug] [-C capQcoef] [-r reg] [-f in.fa] [-l list] [-M capMapQ] [-Q minBaseQ] [-q minMapQ] in.bam [in2.bam [...]]$ samtools mpileup -f genome.fasta abc.bam >abc.txt

$ samtools mpileup -gSDf genome.fasta abc.bam >abc.bcf

$ samtools mpileup -guSDf genome.fasta abc.bam | \

          bcftools view -cvNg - >abc.vcf

mpileup不使用-u或-g參數時(shí)，則不生成二進(jìn)制的bcf文件，而生成一個(gè)文本文件(輸出到標準輸出)。該文本文件統計了參考序列中每個(gè)堿基位點(diǎn)的比對情況；該文件每一行代表了參考序列中某一個(gè)堿基位點(diǎn)的比對結果。比如：

scaffold_1      2841    A      11      ,,,...,....    BHIGDGIJ?FF

scaffold_1      2842    C      12      ,$,,...,....^I. CFGEGEGGCFF+

scaffold_1      2843    G      11      ,,...,.....    FDDDDCD?DD+

scaffold_1      2844    G      11      ,,...,.....    FA?AAAA

scaffold_1      2845    G      11      ,,...,.....    F656666166*

scaffold_1      2846    A      11      ,,...,.....    (1.1111)11*

scaffold_1      2847    A      11      ,,+9acggtgaag.+9ACGGTGAAT.+9ACGGTGAAG.+9ACGGTGAAG,+9acggtgaag.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG      %.+....-..)

scaffold_1      2848    N      11      agGGGgGGGGG    !!$!!!!!!!!

scaffold_1      2849    A      11      c$,...,.....    !0000000000

scaffold_1      2850    A      10      ,...,.....      353333333

mpileup生成的結果包含6行：參考序列名；位置；參考堿基；比對上的reads數；比對情況；比對上的堿基的質(zhì)量。其中第5列比較復雜,解釋如下：

1 ‘.’代表與參考序列正鏈匹配。

2 ‘,’代表與參考序列負鏈匹配。

3 ‘ATCGN’代表在正鏈上的不匹配。

4 ‘a(chǎn)tcgn’代表在負鏈上的不匹配。

5 ‘*’代表模糊堿基

6 ‘^’代表匹配的堿基是一個(gè)read的開(kāi)始；’^'后面緊跟的ascii碼減去33代表比對質(zhì)量；這兩個(gè)符號修飾的是后面的堿基，其后緊跟的堿基(.,ATCGatcgNn)代表該read的第一個(gè)堿基。

7 ‘$’代表一個(gè)read的結束，該符號修飾的是其前面的堿基。

8 正則式’\+[0-9]+[ACGTNacgtn]+’代表在該位點(diǎn)后**入的堿基；比如上例中在scaffold_1的2847后**入了9個(gè)長(cháng)度的堿基acggtgaag。表明此處極可能是indel。

9 正則式’-[0-9]+[ACGTNacgtn]+’代表在該位點(diǎn)后缺失的堿基；

12. 使用bcftools

bcftools和samtools類(lèi)似，用于處理vcf(variant call format)文件和bcf(binary call format)文件。前者為文本文件，后者為其二進(jìn)制文件。

bcftools使用簡(jiǎn)單，最主要的命令是view命令，其次還有index和cat等命令。index和cat命令和samtools中類(lèi)似。此處主講使用view命令來(lái)進(jìn)行SNP和Indel calling。該命令的使用方法和例子為：

$ bcftools view [-AbFGNQSucgv] [-D seqDict] [-l listLoci] [-s listSample]

          [-i gapSNPratio] [-t mutRate] [-p varThres] [-P prior]

          [-1 nGroup1] [-d minFrac] [-U nPerm] [-X permThres]

          [-T trioType] in.bcf [region]

$ bcftools view -cvNg abc.bcf >snp_indel.vcf

生成的結果文件為vcf格式，有10列，分別是：1 參考序列名；2 varianti所在的left-most位置；3 variant的ID（默認未設置，用’.'表示)；4 參考序列的allele；5 variant的allele(有多個(gè)alleles，則用’,'分隔);6 variant/reference QUALity;7 FILTers applied;8 variant的信息，使用分號隔開(kāi)；9 FORMAT of the genotype fields, separated by colon (optional)； 10 SAMPLE genotypes and per-sample information (optional)。

例如：

scaffold_1      2847    .      A      AACGGTGAAG      194    .      INDEL;DP=11;VDB=0.0401;AF1=1;AC1=2;DP4=0,0,8,3;MQ=35;FQ=-67.5  GT:PL:GQ        1/1:235,33,0:63

scaffold_1      3908    .      G      A      111    .      DP=13;VDB=0.0085;AF1=1;AC1=2;DP4=0,0,5,7;MQ=42;FQ=-63  GT:PL:GQ        1/1:144,36,0:69

scaffold_1      4500    .      A      G      31.5    .      DP=8;VDB=0.0034;AF1=1;AC1=2;DP4=0,0,1,3;MQ=42;FQ=-39    GT:PL:GQ        1/1:64,12,0:21

scaffold_1      4581    .      TGGNGG  TGG    145    .      INDEL;DP=8;VDB=0.0308;AF1=1;AC1=2;DP4=0,0,0,8;MQ=42;FQ=-58.5    GT:PL:GQ        1/1:186,24,0:45

scaffold_1      4644    .      G      A      195    .      DP=21;VDB=0.0198;AF1=1;AC1=2;DP4=0,0,10,10;MQ=42;FQ=-87 GT:PL:GQ        1/1:228,60,0:99

scaffold_1      4827    .      NACAAAGA        NA      4.42    .      INDEL;DP=1;AF1=1;AC1=2;DP4=0,0,1,0;MQ=40;FQ=-37.5      GT:PL:GQ        0/1:40,3,0:3

scaffold_1      4854    .      A      G      48      .      DP=6;VDB=0.0085;AF1=1;AC1=2;DP4=0,0,2,1;MQ=41;FQ=-36    GT:PL:GQ        1/1:80,9,0:16

scaffold_1      5120    .      A      G      85      .      DP=8;VDB=0.0355;AF1=1;AC1=2;DP4=0,0,5,3;MQ=42;FQ=-51    GT:PL:GQ        1/1:118,24,0:45

第8列中顯示了對variants的信息描述，比較重要，其中的 Tag 的描述如下：

Tag Format Description

AF1 double Max-likelihood estimate of the site allele frequency (AF) of the first ALT allele

DP int Raw read depth (without quality filtering)

DP4 int[4] # high-quality reference forward bases, ref reverse, alternate for and alt rev bases

FQ int Consensus quality. Positive: sample genotypes different; negative: otherwise

MQ int Root-Mean-Square mapping quality of covering reads

PC2 int[2] Phred probability of AF in group1 samples being larger (,smaller) than in group2

PCHI2 double Posterior weighted chi^2 P-value between group1 and group2 samples

PV4 double[4] P-value for strand bias, baseQ bias, mapQ bias and tail distance bias

QCHI2 int Phred-scaled PCHI2

RP int # permutations yielding a smaller PCHI2

CLR int Phred log ratio of genotype likelihoods with and without the trio/pair constraint

UGT string Most probable genotype configuration without the trio constraint

CGT string Most probable configuration with the trio constraint

使用bcftools得到variant calling結果后。需要對結果再次進(jìn)行過(guò)濾。主要依據比對結果中第8列信息。其中的 DP4 一行尤為重要，提供了4個(gè)數據：1 比對結果和正鏈一致的reads數、2 比對結果和負鏈一致的reads數、3 比對結果在正鏈的variant上的reads數、4 比對結果在負鏈的variant上的reads數?？梢栽O定 （value3 + value4）大于某一閾值，才算是variant。比如：

$ perl -ne 'print $_ if /DP4=(\d+),(\d+),(\d+),(\d+)/ && ($3+$4)>=10 && ($3+$4)/($1+$2+$3+$4)>=0.8' snp_indel.vcf >snp_indel.final.vcf