上次我們整理到bwa比對(duì)后得到bam文件，下一步我們要通過(guò)GATK流程從bam文件中call variant。
一、使用GATK前須知事項(xiàng)：
對(duì)GATK的測(cè)試主要使用的是人類(lèi)全基因組和外顯子組的測(cè)序數(shù)據(jù)，而且全部是基于illumina數(shù)據(jù)格式，目前還沒(méi)有提供其他格式文件（如Ion Torrent）或者實(shí)驗(yàn)設(shè)計(jì)（RNA-Seq）的分析方法。
GATK是一個(gè)應(yīng)用于前沿科學(xué)研究的軟件，不斷在更新和修正，因此，在使用GATK進(jìn)行變異檢測(cè)時(shí)，最好是下載最新的版本。
在GATK使用過(guò)程中，有些步驟需要用到已知變異信息，對(duì)于這些已知變異，GATK只提供了人類(lèi)的已知變異信息 ，可以在GATK的FTP站點(diǎn)下載（GATK resource bundle）。如果要研究的不是人類(lèi)基因組，需要自行構(gòu)建已知變異，GATK提供了詳細(xì)的構(gòu)建方法。
GATK在進(jìn)行BQSR和VQSR的過(guò)程中會(huì)使用到R軟件繪制一些圖，因此，在運(yùn)行GATK之前最好先檢查一下是否正確安裝了R和所需要的包，所需要的包大概包括ggplot2、gplots、bitops、caTools、colorspace、gdata、gsalib、reshape、RColorBrewer等。如果畫(huà)圖時(shí)出現(xiàn)錯(cuò)誤，會(huì)提示需要安裝的包的名稱(chēng)。
第一步：原始數(shù)據(jù)的處理：對(duì)原始下機(jī)fastq文件進(jìn)行過(guò)濾和比對(duì)（mapping）對(duì)于Illumina下機(jī)數(shù)據(jù)推薦使用bwa進(jìn)行mapping。前邊已講。
bwa中 -r str：定義頭文件?！瓳RG\tID:foo\tSM:bar’，如果在此步驟不進(jìn)行頭文件定義，在后續(xù)GATK分析中還是需要重新增加頭文件。GATK2.0以上版本將不再支持無(wú)頭文件的變異檢測(cè)。
ID str：輸入reads集ID號(hào)；LB：read集文庫(kù)名；PL：測(cè)序平臺(tái)（illunima或solid）；PU：測(cè)序平臺(tái)下級(jí)單位名稱(chēng)（run的名稱(chēng)）；SM：樣本名稱(chēng)
對(duì)于最后得到的sam文件，將比對(duì)上的結(jié)果提取出來(lái)（awk即可處理），即可直接用于GATK的分析。
注意：由于GATK在下游的snpcalling時(shí)，是按染色體進(jìn)行callsnp的。因此，在準(zhǔn)備原始sam文件時(shí)，可以先按染色體將文件分開(kāi)，這樣會(huì)提高運(yùn)行速度。但是當(dāng)數(shù)據(jù)量不足時(shí)，可能會(huì)影響后續(xù)的VQSR分析，這是需要注意的。
對(duì)sam文件進(jìn)行進(jìn)行重新排序（reorder）
由BWA生成的sam文件時(shí)按字典式排序法進(jìn)行的排序（lexicographically）進(jìn)行排序的（chr10，chr11…chr19，chr1，chr20…chr22，chr2，chr3…chrM，chrX，chrY），但是GATK在進(jìn)行callsnp的時(shí)候是按照染色體組型（karyotypic）進(jìn)行的（chrM，chr1，chr2…chr22，chrX，chrY），因此要對(duì)原始sam文件進(jìn)行reorder?？梢允褂胮icard-tools中的ReorderSam完成。
對(duì)bam文件進(jìn)行sort排序處理 一步是將sam文件中同一染色體對(duì)應(yīng)的條目按照坐標(biāo)順序從小到大進(jìn)行排序?？梢允褂胮icard-tools中SortSam完成。
第一步： Duplicates Marking
在制備文庫(kù)的過(guò)程中，由于PCR擴(kuò)增過(guò)程中會(huì)存在一些偏差，也就是說(shuō)有的序列會(huì)被過(guò)量擴(kuò)增。這樣，在比對(duì)的時(shí)候，這些過(guò)量擴(kuò)增出來(lái)的完全相同的序列就會(huì)比對(duì)到基因組的相同位置。而這些過(guò)量擴(kuò)增的reads并不是基因組自身固有序列，不能作為變異檢測(cè)的證據(jù)，因此，要盡量去除這些由PCR擴(kuò)增所形成的duplicates，這一步可以使用picard-tools來(lái)完成。
去重復(fù)的過(guò)程是給這些序列設(shè)置一個(gè)flag以標(biāo)志它們，方便GATK的識(shí)別。還可以設(shè)置 REMOVE_DUPLICATES=true 來(lái)丟棄duplicated序列。對(duì)于是否選擇標(biāo)記或者刪除，對(duì)結(jié)果應(yīng)該沒(méi)有什么影響，GATK官方流程里面給出的例子是僅做標(biāo)記不刪除。這里定義的重復(fù)序列是這樣的：如果兩條reads具有相同的長(zhǎng)度而且比對(duì)到了基因組的同一位置，那么就認(rèn)為這樣的reads是由PCR擴(kuò)增而來(lái)，就會(huì)被GATK標(biāo)記。
最主要目的就是盡量減小文庫(kù)構(gòu)建時(shí)引入文庫(kù)的PCR bias
標(biāo)記后對(duì)結(jié)果文件生成索引文件。
GATK=/home/jmzeng/biosoft/gatk4/gatk-4.0.6.0/gatkref=/public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fastasnp=/public/biosoft/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gzindel=/public/biosoft/GATK/resources/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gzforsamplein{7E5239.L1,7E5240,7E5241.L1}doecho $sample$GATK --java-options"-Xmx20G -Djava.io.tmpdir=./"MarkDuplicates \? ? -I $sample.bam \? ? -O ${sample}_marked.bam \? ? -M $sample.metrics \1>${sample}_log.mark2>&1samtools index ${sample}_marked_fixed.bam
第二步：Local realignment around indels
這一步的目的就是將比對(duì)到indel附近的reads進(jìn)行局部重新比對(duì)，將比對(duì)的錯(cuò)誤率降到最低。一般來(lái)說(shuō)，絕大部分需要進(jìn)行重新比對(duì)的基因組區(qū)域，都是因?yàn)椴迦?缺失的存在，因?yàn)樵趇ndel附近的比對(duì)會(huì)出現(xiàn)大量的堿基錯(cuò)配，這些堿基的錯(cuò)配很容易被誤認(rèn)為SNP。還有，在比對(duì)過(guò)程中，比對(duì)算法對(duì)于每一條read的處理都是獨(dú)立的，不可能同時(shí)把多條reads與參考基因組比對(duì)來(lái)排錯(cuò)。因此，即使有一些reads能夠正確的比對(duì)到indel，但那些恰恰比對(duì)到indel開(kāi)始或者結(jié)束位置的read也會(huì)有很高的比對(duì)錯(cuò)誤率，這都是需要重新比對(duì)的。
主要分為兩步：
通過(guò)運(yùn)行RealignerTargetCreator來(lái)確定要進(jìn)行重新比對(duì)的區(qū)域。
java-jarGenomeAnalysisTK.jar-Rhg19.fa-TRealignerTargetCreator-Ihg19.reorder.sort.addhead.dedup_04.bam-ohg19.dedup.realn_06.intervals-knownMills_and_1000G_gold_standard.indels.hg38.vcf-known1000G_phase1.indels.hg38.vcf參數(shù)說(shuō)明：-R： 參考基因組；-T： 選擇的GATK工具；-I：? 輸入上一步所得bam文件；-o： 輸出的需要重新比對(duì)的基因組區(qū)域結(jié)果；-maxInterval:允許進(jìn)行重新比對(duì)的基因組區(qū)域的最大值，不能太大，太大耗費(fèi)會(huì)很長(zhǎng)時(shí)間，默認(rèn)值500；-known:已知的可靠的indel位點(diǎn)，指定已知的可靠的indel位點(diǎn)，重比對(duì)將主要圍繞這些位點(diǎn)進(jìn)行，對(duì)于人類(lèi)基因組數(shù)據(jù)而言，可以直接指定GATKresourcebundle里面的indel文件（必須是vcf文件）。
通過(guò)運(yùn)行IndelRealigner在這些區(qū)域內(nèi)進(jìn)行重新比對(duì)
-R hg19.fa-T IndelRealigner-targetIntervals hg19.dedup.realn_06.intervals-I hg19.reorder.sort.addhead.dedup_04.bam-o hg19.dedup.realn_07.bam-known Mills_and_1000G_gold_standard.indels.hg38.vcf-known1000G_phase1.indels.hg38.vcf
注意：1. 第一步和第二步中使用的輸入文件（bam文件）、參考基因組和已知indel文件必須是相同的文件。
2. 當(dāng)在相同的基因組區(qū)域發(fā)現(xiàn)多個(gè)indel存在時(shí)，這個(gè)工具會(huì)從其中選擇一個(gè)最有可能存在比對(duì)錯(cuò)誤的indel進(jìn)行重新比對(duì)，剩余的其他indel不予考慮。
第三步：.Base quality score recalibration
這一步是對(duì)bam文件里reads的堿基質(zhì)量值進(jìn)行重新校正，使最后輸出的bam文件中reads中堿基的質(zhì)量值能夠更加接近真實(shí)的與參考基因組之間錯(cuò)配的概率。這一步適用于多種數(shù)據(jù)類(lèi)型，包括illunima、solid、454、CG等數(shù)據(jù)格式。在GATK2.0以上版本中還可以對(duì)indel的質(zhì)量值進(jìn)行校正，這一步對(duì)indel calling非常有幫助。
舉例說(shuō)明，在reads堿基質(zhì)量值被校正之前，我們要保留質(zhì)量值在Q25以上的堿基，但是實(shí)際上質(zhì)量值在Q25的這些堿基的錯(cuò)誤率在1%，也就是說(shuō)質(zhì)量值只有Q20，這樣就會(huì)對(duì)后續(xù)的變異檢測(cè)的可信度造成影響。還有，在邊合成邊測(cè)序的測(cè)序過(guò)程中，在reads末端堿基的錯(cuò)誤率往往要比起始部位更高。另外，AC的質(zhì)量值往往要低于TG。BQSR的就是要對(duì)這些質(zhì)量值進(jìn)行校正。
利用工具BaseRecalibrator，根據(jù)一些known sites，生成一個(gè)校正質(zhì)量值所需要的數(shù)據(jù)文件，GATK網(wǎng)站以“.grp”為后綴命名。
-R$ref\? ? -I${sample}_marked_fixed.bam? \? ? --known-sites$snp\? ? --known-sites$indel\? ? -O${sample}_recal.table \? ? 1>${sample}_log.recal 2>&1
第四步：ApplyBQSR
利用已知snp數(shù)據(jù)庫(kù)，通過(guò)機(jī)器學(xué)習(xí)方法調(diào)整堿基質(zhì)量分?jǐn)?shù)
-R$ref\? ? -I${sample}_marked_fixed.bam? \? ? -bqsr${sample}_recal.table \? ? -O${sample}_bqsr.bam \? ? 1>${sample}_log.ApplyBQSR? 2>&1
這時(shí)候，GATK流程文件準(zhǔn)備工作結(jié)束，完成了variants calling所需要的所有準(zhǔn)備，生成了用于下一步變異檢測(cè)的bam文件。
第五步：Variant Calling
GATK在這一步里面提供了兩個(gè)工具進(jìn)行變異檢測(cè)——UnifiedGenotyper和HaplotypeCaller
HaplotypeCaller不支持Reduce之后的bam文件，因此，當(dāng)選擇使用HaplotypeCaller進(jìn)行變異檢測(cè)時(shí)，不需要進(jìn)行Reduce reads
-R$ref\? ? -I${sample}_bqsr.bam \? ? ? --dbsnp$snp\? ? ? -O${sample}_raw.vcf \? ? ? 1>${sample}_log.HC 2>&1
對(duì)輸入的bam文件中的所有樣本進(jìn)行變異檢測(cè)，最后生成一個(gè)vcf文件，vcf文件中會(huì)包含所有樣本的變異位點(diǎn)和基因型信息。但是現(xiàn)在所得到的結(jié)果是最原始的、沒(méi)有經(jīng)過(guò)任何過(guò)濾和校正的Variants集合。這一步產(chǎn)生的變異位點(diǎn)會(huì)有很高的假陽(yáng)性，尤其是indel，因此，必須要進(jìn)行進(jìn)一步的篩選過(guò)濾。這一步還可以指定對(duì)基因組的某一區(qū)域進(jìn)行變異檢測(cè)，只需要增加一個(gè)參數(shù) -L：target_interval.list，target_interval.list 格式是bed格式文件。
bed文件格式，三列必須
TCGAGA#對(duì)應(yīng)bed文件的坐標(biāo)應(yīng)為#chrome start endchr1? ? ? ? ? ? 0? ? 5
注意：GATK進(jìn)行變異檢測(cè)的時(shí)候，是按照染色體排序順序進(jìn)行的（先call chr1，然后chr2，然后chr3…最后chrY），并非多條染色體并行檢測(cè)的，因此，如果數(shù)據(jù)量比較大的話，建議分染色體分別進(jìn)行，對(duì)性染色體的變異檢測(cè)可以同常染色體方法。
大多數(shù)參數(shù)的默認(rèn)值可以滿足大多數(shù)研究的需求，因此，在做變異檢測(cè)過(guò)程中，如果對(duì)參數(shù)意義不是很明確，不建議修改。
第六步： 對(duì)原始變異檢測(cè)結(jié)果進(jìn)行過(guò)濾（hard filter and VQSR）
這一步的目的就是對(duì)上一步call出來(lái)的變異位點(diǎn)進(jìn)行過(guò)濾，去掉不可信的位點(diǎn)。這一步可以有兩種方法，一種是通過(guò)GATK的VariantFiltration，另一種是通過(guò)GATK的VQSR（變異位點(diǎn)質(zhì)量值重新校正）進(jìn)行過(guò)濾。
GATK是推薦使用VQSR的，但使用VQSR數(shù)據(jù)量一定要達(dá)到要求，數(shù)據(jù)量太小無(wú)法使用高斯模型。在使用VAQR時(shí)，indel和snp要分別進(jìn)行。
VQSR原理介紹：
這個(gè)模型是根據(jù)已有的真實(shí)變異位點(diǎn)（人類(lèi)基因組一般使用HapMap3中的位點(diǎn)，以及這些位點(diǎn)在Omni 2.5M SNP芯片中出現(xiàn)的多態(tài)位點(diǎn)）來(lái)訓(xùn)練，最后得到一個(gè)訓(xùn)練好的能夠很好的評(píng)估真?zhèn)蔚腻e(cuò)誤評(píng)估模型，可以叫他適應(yīng)性錯(cuò)誤評(píng)估模型。這個(gè)適應(yīng)性的錯(cuò)誤評(píng)估模型可以應(yīng)用到call出來(lái)的原始變異集合中已知的變異位點(diǎn)和新發(fā)現(xiàn)的變異位點(diǎn)，進(jìn)而去評(píng)估每一個(gè)變異位點(diǎn)發(fā)生錯(cuò)誤的概率，最終會(huì)給出一個(gè)得分。這個(gè)得分最后會(huì)被寫(xiě)入vcf文件的INFO信息里，叫做VQSLOD，就是在訓(xùn)練好的混合高斯模型下，一個(gè)位點(diǎn)是真實(shí)的概率比上這個(gè)位點(diǎn)可能是假陽(yáng)性的概率的log odds ratio（對(duì)數(shù)差異比），因此，可以定性的認(rèn)為，這個(gè)值越大就越好。
VQSR主要分兩個(gè)步驟，這兩個(gè)步驟會(huì)使用兩個(gè)不同的工具：VariantRecalibrator和ApplyRecalibration。
VariantRecalibrator：通過(guò)大量的高質(zhì)量的已知變異集合的各個(gè)注釋?zhuān)òê芏喾N，后面介紹）的值來(lái)創(chuàng)建一個(gè)高斯混合模型，然后用于評(píng)估所有的變異位點(diǎn)。這個(gè)文件最后將生成一個(gè)recalibration文件。
原理簡(jiǎn)單介紹： 這個(gè)模型首先要拿到真實(shí)變異數(shù)據(jù)集和上一步驟中得到的原始變異數(shù)據(jù)集的交集，然后對(duì)這些SNP值相對(duì)于具體注釋信息的分布情況進(jìn)行模擬，將這些變異位點(diǎn)進(jìn)行聚類(lèi)，最后根據(jù)聚類(lèi)結(jié)果賦予所有變異位點(diǎn)相應(yīng)的VQSLOD值。越接近聚類(lèi)核心的變異位點(diǎn)得到的VQSLOD值越高。
ApplyRecalibration：這一步將模型的各個(gè)參數(shù)應(yīng)用于原始vcf文件中的每一個(gè)變異位點(diǎn)，這時(shí)，每一個(gè)變異位點(diǎn)的注釋信息列中都會(huì)出現(xiàn)一個(gè)VQSLOD值，然后模型會(huì)根據(jù)這個(gè)值對(duì)變異位點(diǎn)進(jìn)行過(guò)濾，過(guò)濾后的信息會(huì)寫(xiě)在vcf文件的filter一列中。
原理簡(jiǎn)單介紹：在VariantRecalibrator這一步中，每個(gè)變異位點(diǎn)已經(jīng)得到了一個(gè)VQSLOD值了，同時(shí)，這些LOD值在訓(xùn)練集里也進(jìn)行了排序。當(dāng)你在這一步中設(shè)置一個(gè)tranche sensitivity 的閾值（這個(gè)閾值一般是一個(gè)百分?jǐn)?shù)，如設(shè)置成99%），那么，如果LOD值從大到小排序的話，這個(gè)程序就會(huì)認(rèn)為在這個(gè)訓(xùn)練集中，LOD值在前99%的是可信的，當(dāng)這個(gè)值低于這個(gè)閾值，就認(rèn)為是錯(cuò)誤的。最后，程序就會(huì)用這個(gè)標(biāo)準(zhǔn)來(lái)過(guò)濾上一步call出來(lái)的原始變異集合。如果LOD值超過(guò)這個(gè)閾值，在filter那一列就會(huì)顯示PASS，如果低于這個(gè)值就會(huì)被過(guò)濾掉，但是這些位點(diǎn)仍然會(huì)顯示在結(jié)果里面，只不過(guò)會(huì)在filter那一列標(biāo)示出他所屬于的tranche sensitivity 的名稱(chēng)。在設(shè)置tranche sensitivity 的閾值時(shí)，要兼顧敏感度和質(zhì)量值。
tranche值的設(shè)定
前面提到了，這個(gè)值得設(shè)定是用來(lái)在后續(xù)的ApplyRecalibration中如何根據(jù)這個(gè)閾值來(lái)過(guò)濾變異位點(diǎn)的，也就是說(shuō)，如果這個(gè)值設(shè)定的比較高的話，那么最后留下來(lái)的變異位點(diǎn)就會(huì)多，但同時(shí)假陽(yáng)性的位點(diǎn)也會(huì)相應(yīng)增加；如果設(shè)定的低的話，雖然假陽(yáng)性會(huì)減少，但是會(huì)丟失很多真實(shí)的位點(diǎn)。因此，跟選擇注釋時(shí)一樣，可以run兩遍VariantRecalibrator，第一遍的時(shí)候多寫(xiě)幾個(gè)閾值，第一遍跑完之后看結(jié)果，看那個(gè)閾值好，選擇一個(gè)最好的閾值，再run一遍VariantRecalibrator。至于說(shuō)怎么區(qū)分好壞，有幾個(gè)標(biāo)準(zhǔn)：
看結(jié)果中已知變異位點(diǎn)與新發(fā)現(xiàn)變異位點(diǎn)之間的比例，這個(gè)比例不要太大，因?yàn)榇蠖鄶?shù)新發(fā)現(xiàn)的變異都是假陽(yáng)性，如果太多的話，可能假陽(yáng)性的比例就比較大；
看保留的變異數(shù)目，這個(gè)就要根據(jù)具體的需求進(jìn)行選擇了。
看TI/TV值，對(duì)于人類(lèi)全基因組，這個(gè)值應(yīng)該在2.15左右，對(duì)于外顯子組，這個(gè)值應(yīng)該在3.2左右，不要太小或太大，越接近這個(gè)數(shù)值越好，這個(gè)值如果太小，說(shuō)明可能存在比較多的假陽(yáng)性。
千人中所選擇的tranche值是99
注意：Indel不支持tranche值的選擇，另外，一部分注釋類(lèi)型在做indel的校正時(shí)也不支持，具體信息可以詳查GATK網(wǎng)站。
當(dāng)數(shù)據(jù)量太小時(shí)，可能高斯模型不會(huì)運(yùn)行，因?yàn)樽儺愇稽c(diǎn)數(shù)滿足不了模型的統(tǒng)計(jì)需求。這時(shí)候可以通過(guò)降低--maxGaussian的值，讓程序運(yùn)行。這個(gè)值表示的是程序?qū)⒆儺愇稽c(diǎn)分成的最大的組數(shù)，降低這個(gè)值讓程序把變異位點(diǎn)聚類(lèi)到更少的組里面，使每個(gè)組中的變異位點(diǎn)數(shù)增加來(lái)滿足統(tǒng)計(jì)需求，但是這樣做降低程序分辨真?zhèn)蔚哪芰ΑＲ虼?，在運(yùn)行程序的時(shí)候，要對(duì)各方面進(jìn)行權(quán)衡。
過(guò)濾后生成的vcf文件的格式說(shuō)明，即每一列所代表的的內(nèi)容，可參考下面的網(wǎng)站，有詳細(xì)的說(shuō)明
http://www.broadinstitute.org/gatk/guide/article?id=1268
第七步：初步注釋分析利用GATK中的VariantEval來(lái)完成，或者用annovar注釋
ANNOVAR是由王凱編寫(xiě)的一個(gè)注釋軟件，可以對(duì)SNP和indel進(jìn)行注釋?zhuān)部梢赃M(jìn)行變異的過(guò)濾篩選。
ANNOVAR能夠利用最新的數(shù)據(jù)來(lái)分析各種基因組中 的遺傳變異。主要包含三種不同的注釋方法，Gene-based Annotation（基于基因的注釋?zhuān)?、Region-based Annotation（基于區(qū)域的注釋?zhuān)ilter-based Annotation（基于篩選的注釋?zhuān)?。ANNOVAR由Perl編寫(xiě)。
優(yōu)點(diǎn)：提供多個(gè)數(shù)據(jù)可直接下載、支持多種格式、注釋直觀；缺點(diǎn)：沒(méi)有數(shù)據(jù)庫(kù)的物種無(wú)法注釋。
│? annotate_variation.pl#主程序，功能包括下載數(shù)據(jù)庫(kù)，三種不同的注釋│? coding_change.pl#可用來(lái)推斷蛋白質(zhì)序列│? convert2annovar.pl#將多種格式轉(zhuǎn)為.avinput的程序│? retrieve_seq_from_fasta.pl#用于自行建立其他物種的轉(zhuǎn)錄本│? table_annovar.pl#注釋程序，可一次性完成三種類(lèi)型的注釋│? variants_reduction.pl#可用來(lái)更靈活地定制過(guò)濾注釋流程│├─example#存放示例文件│└─humandb#人類(lèi)注釋數(shù)據(jù)庫(kù)
ANNOVAR下載數(shù)據(jù)庫(kù)
命令示例
# -buildver 表示version# -downdb 下載數(shù)據(jù)庫(kù)的指令# -webfrom annovar 從annovar提供的鏡像下載，不加此參數(shù)將尋找數(shù)據(jù)庫(kù)本身的源# humandb/ 存放于humandb/目錄下
輸出的csv文件將包含輸入的5列主要信息以及各個(gè)數(shù)據(jù)庫(kù)里的注釋?zhuān)送猓瑃able_annoval.pl可以直接對(duì)vcf文件進(jìn)行注釋?zhuān)ú恍枰D(zhuǎn)換格式），注釋的內(nèi)容將會(huì)放在vcf文件的“INFO”那一欄。
說(shuō)明：本文主要說(shuō)明wes分析流程的理論，代碼僅作示例。
參考文獻(xiàn)
1.http://www.itdecent.cn/p/49d035b121b8
2.https://blog.csdn.net/zhu_si_tao/article/details/53321374
3.https://www.plob.org/article/9976.html
作者：鳳凰_0949
鏈接：http://www.itdecent.cn/p/f9ac13d3c938
來(lái)源：簡(jiǎn)書(shū)
簡(jiǎn)書(shū)著作權(quán)歸作者所有，任何形式的轉(zhuǎn)載都請(qǐng)聯(lián)系作者獲得授權(quán)并注明出處。
