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

PCHI英文含義

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

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

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

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

VCF文件參數(shù)解讀

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

    例如：以##FILTER開(kāi)頭的行，表示注釋VCF文件當(dāng)中第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     記錄染色體編號(hào)

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

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

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

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

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

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

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

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

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

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

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

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

    #QD：通過(guò)深度來(lái)評(píng)估一個(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的數(shù)目；AF(Allele Frequency) 表示Allele的頻率； AN(Allele Number) 表示Allele的總數(shù)目。

     對(duì)于1個(gè)diploid sample（雙倍體）而言：則基因型 0/1 表示sample為雜合子，Allele數(shù)為1 (雙倍體的sample在該位點(diǎn)只有1個(gè)等位基因發(fā)生了突變)，Allele的頻率為0.5 (雙倍體的sample在該位點(diǎn)只有50%的等位基因發(fā)生了突變)，總的Allele為2；                                                                    基因型 1/1 則表示sample為純合的，Allele數(shù)為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ì)量，負(fù)值表示支持變異的堿基質(zhì)量值不及支持參考基因組的，

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

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

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

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

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

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

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

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

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

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

    #HaplotypeScore    分?jǐn)?shù)越高越可能出現(xiàn)錯(cuò)誤。Higher scores are indicative of regions with bad alignments, typically leading to artifactual SNP and indel calls.

    #SOR：也是一個(gè)用來(lái)評(píng)估是否存在鏈偏向性的參數(shù)，相當(dāng)于FS的升級(jí)版。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允許的最大數(shù)量

    #G3：ML 評(píng)估基因型出現(xiàn)的頻率

    #HWE：chi^2基于HWE的測(cè)試p值和G3

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

    #UGT：最可能不受限制的三種基因型結(jié)構(gòu)

    #CGT：最可能受限制三種基因型的結(jié)構(gòu)

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

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

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

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

    #QCHI2：Phred scaled PCHI2

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

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

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

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

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

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

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

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

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

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

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

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

     PL ： 指定的三種基因型的似然值。這三種指定的基因型為(0/0,0/1,1/1)，這三種基因型的概率總和為1。數(shù)值越小代表基因型越可靠，最小的數(shù)字對(duì)應(yīng)的基因型判讀為該樣品的最可能的基因型。比如最后一列285,0,105，分別對(duì)應(yīng)基因型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文件（通常是短序列比對(duì)工具如bwa,bowtie2,hisat2,tophat2等等產(chǎn)生的，具體格式可以在消息框輸入“SAM”查看）的工具合集，包含有許多命令。以下是常用命令的介紹。

1.View

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

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

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

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

默認(rèn)情況下不加 region，則是輸出所有的 region.

options:

-b output BAM

  默認(rèn)下輸出是 SAM 格式文件，該參數(shù)設(shè)置輸出 BAM 格式

-h print header for the SAM output

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

-H print header only (no alignments)

  僅僅輸出文件的頭

-S input is SAM

  默認(rèn)下輸入是 BAM 文件，若是輸入是 SAM 文件，則最好加該參數(shù)，否則有時(shí)候會(huì)報(bào)錯(cuò)。

-u uncompressed BAM output (force -b)

  該參數(shù)的使用需要有-b參數(shù)，能節(jié)約時(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ò)濾功統(tǒng)計(jì)功能

-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]

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

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

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

-q INT  minimum mapping quality [0]

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

例子：

將sam文件轉(zhuǎn)換成bam文件

samtools view -bS abc.sam >abc.bam

BAM轉(zhuǎn)換為SAM

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

提取比對(duì)到參考序列上的比對(duì)結(jié)果

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

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

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

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

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

提取bam文件中比對(duì)到caffold1上的比對(duì)結(jié)果，并保存到sam文件格式

samtools view abc.bam scaffold1 >scaffold1.sam

提取scaffold1上能比對(duì)到30k到100k區(qū)域的比對(duì)結(jié)果

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

根據(jù)fasta文件，將 header 加入到 sam 或 bam 文件中

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

2.Sort

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

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

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

-n 設(shè)定排序方式按short reads的ID排序。默認(rèn)下是按序列在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

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

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

Usage: samtools index [out.index]

例子：

以下兩種命令結(jié)果一樣

$ samtools index abc.sort.bam

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

5.faidx

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

Usage: samtools faidx [ [...]]

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

例子：$ samtools faidx genome.fasta

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

$ samtools faidx genome.fasta scffold_10 >scaffold_10.fasta

6.tview

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

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

Usage: samtools tview [ref.fasta]

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

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

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

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

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

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

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

20～30**；10～20綠色；0～10藍(lán)色。

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

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

7.flagstat

給出BAM文件的比對(duì)結(jié)果

Usage: samtools flagstat

$ samtools flagstat example.bam

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

#總共的reads數(shù)

0 + 0 duplicates

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

#總體上reads的匹配率

11945742 + 0 paired in sequencing

#有多少reads是屬于paired reads

5972871 + 0 read1

#reads1中的reads數(shù)

5972871 + 0 read2

#reads2中的reads數(shù)

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

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

6899708 + 0 with itself and mate mapped

#paired reads中兩條都比對(duì)到參考序列上的reads數(shù)

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

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

469868 + 0 with mate mapped to a different chr

#paired reads中兩條分別比對(duì)到兩條不同的參考序列的reads數(shù)

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

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

8.depth

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

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

-r 后面跟染色體號(hào)（region）

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

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

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

例子

samtools depth accept.bam >depth

9.其他命令

reheader：替換bam文件的頭

$ samtools reheader

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

$ samtools idxstats

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

Usage:  samtools rmdup [-sS] 

-s 對(duì)single-end reads。默認(rèn)情況下，只對(duì)paired-end reads

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

10. 將bam文件轉(zhuǎn)換為fastq文件

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

該網(wǎng)站提供了一個(gè)bam轉(zhuǎn)換為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的安裝文件夾中可以找到。

最常用的參數(shù)有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ù)時(shí)，則不生成二進(jìn)制的bcf文件，而生成一個(gè)文本文件(輸出到標(biāo)準(zhǔn)輸出)。該文本文件統(tǒng)計(jì)了參考序列中每個(gè)堿基位點(diǎn)的比對(duì)情況；該文件每一行代表了參考序列中某一個(gè)堿基位點(diǎn)的比對(duì)結(jié)果。比如：

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生成的結(jié)果包含6行：參考序列名；位置；參考?jí)A基；比對(duì)上的reads數(shù)；比對(duì)情況；比對(duì)上的堿基的質(zhì)量。其中第5列比較復(fù)雜,解釋如下：

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

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

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

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

5 ‘*’代表模糊堿基

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

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

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

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

12. 使用bcftools

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

bcftools使用簡(jiǎn)單，最主要的命令是view命令，其次還有index和cat等命令。index和cat命令和samtools中類似。此處主講使用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

生成的結(jié)果文件為vcf格式，有10列，分別是：1 參考序列名；2 varianti所在的left-most位置；3 variant的ID（默認(rèn)未設(shè)置，用’.'表示)；4 參考序列的allele；5 variant的allele(有多個(gè)alleles，則用’,'分隔);6 variant/reference QUALity;7 FILTers applied;8 variant的信息，使用分號(hào)隔開(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列中顯示了對(duì)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結(jié)果后。需要對(duì)結(jié)果再次進(jìn)行過(guò)濾。主要依據(jù)比對(duì)結(jié)果中第8列信息。其中的 DP4 一行尤為重要，提供了4個(gè)數(shù)據(jù)：1 比對(duì)結(jié)果和正鏈一致的reads數(shù)、2 比對(duì)結(jié)果和負(fù)鏈一致的reads數(shù)、3 比對(duì)結(jié)果在正鏈的variant上的reads數(shù)、4 比對(duì)結(jié)果在負(fù)鏈的variant上的reads數(shù)?？梢栽O(shè)定 （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