• 官網(wǎng)：http://grass.cgs.hku.hk/limx/kggseq/

• kggseq軟件是一款功能強(qiáng)大的測序后下游分析軟件，現(xiàn)簡要列舉其功能如下，以備查詢。

1. 質(zhì)控（Quality Control）

1. 對基因型的質(zhì)控（Genotype QC）
   參數(shù)如下：

   
   

2. 對位點(diǎn)的質(zhì)控（Variant QC）
   參數(shù)如下：

   
   
   image.png
   
   
   image.png
3. 對樣本的質(zhì)控（Sample QC）
   可以使用kggseq軟件把vcf轉(zhuǎn)換為plink格式進(jìn)行質(zhì)控，包括 (pairwise relatedness, mendel error rate and sex check)

java -jar kggseq.jar --vcf-file path-to-file1 --ped-file path/to/file2 <Variants and Genotype QC Settings> --db-merge <Reference genotype set> --o-plink-bed 
plink --bfile kggseq.merged --genome
plink --bfile kggseq --mendel
plink --bfile kggseq --check-sex

• 也可以畫MAF分布圖

java -jar kggseq.jar --vcf-file path-to-file1 --out path/to/prefixname --mafplot
java -jar kggseq.jar --vcf-file path-to-file1 --ped-file path/to/file2  --db-filter hg18_dbsnp138,hg18_dbsnp141 --allele-freq 0,0.01 --o-vcf path/to/file3

2. 過濾（Filtering）

1. 使用--genotype-filter參數(shù)，后面加上數(shù)字，數(shù)字對應(yīng)意義如下：

   
   
   image.png

• 使用--double-hit-gene-trio-filter 或 --double-hit-gene-phased-filter參數(shù)提取復(fù)合雜合或者隱形突變位點(diǎn)

java -jar ./kggseq.jar --buildver hg19 --vcf-file path/to/file1 --ped-file path/to/file2 --db-gene refgene --gene-feature-in 0,1,2,3,4,5,6 --db-filter 1kg201204,dbsnp137,ESP6500AA,ESP6500EA  --rare-allele-freq 0.033 --db-score dbnsfp --mendel-causing-predict best --double-hit-gene-trio-filter (OR --double-hit-gene-phased-filter)

2. 使用基因特征過濾（Gene feature filtering）
   參數(shù)：--db-gene & --gene-feature-in （兩個參數(shù)分別選擇數(shù)據(jù)庫和基因特征）

java -jar ./kggseq.jar  --buildver hg19 --vcf-file path/to/file1 --ped-file path/to/file2 --dgv-cnv-annot --splicing 6 --db-gene refgene,gencode,knowngene --gene-feature-in 0,1,2,3,4,5,6 

基因特征代號如下：

image.png

3. 使用公共數(shù)據(jù)庫的基因型頻率進(jìn)行過濾（--db-filter）

java -jar ./kggseq.jar --buildver hg19 --vcf-file path/to/file1 --ped-file path/to/file2  --db-filter 1kg201204,dbsnp141,ESP6500AA,ESP6500EA  --rare-allele-freq 0.005 --db-filter-hard dbsnp138nf

也可以對本地?cái)?shù)據(jù)進(jìn)行頻率過濾（--local-filter-hard）

4. 過濾掉重復(fù)區(qū)域（--superdup-filter）或者注釋這些區(qū)域（--superdup-annot）

java -jar ./kggseq.jar --buildver hg19 --vcf-file path/to/file1 --ped-file path/to/file2 --superdup-filter

5. 過濾掉indel，某些基因或某些區(qū)域

• --ignore-indel
• --ignore-snv
• --genes-out
• --regions-out

6. 過濾掉假基因等其他因素造成的變異位點(diǎn)過多的基因（--gene-var-filter）

java -jar ./kggseq.jar --buildver hg19 --vcf-file path/to/file1 --ped-file path/to/file2 [options to filter out common and neutral variants] --gene-var-filter 4

7. 過濾掉實(shí)驗(yàn)組/對照組獨(dú)有的位點(diǎn)（--filter-model case-unique or --filter-model control-unique）

java -jar ./kggseq.jar --buildver hg19  --vcf-file path/to/file1 --ped-file path/to/file2 --filter-model case-unique 

3. 注釋（Annotation）

• 能注釋的信息包括：

1. SNP rs號（--rsid）
2. 選擇性剪切位點(diǎn)（--scsnv-annot）
3. 該位點(diǎn)是否位于CNV區(qū)域（--dgv-cnv-annot）
4. 該位點(diǎn)是否與感興趣的基因在同一蛋白通路上（--candi-list gene1,gene2,...,geneN --ppi-annot ppiDatabase ）
5. 該位點(diǎn)是否與感興趣的基因在同一基因通路上（--candi-list gene1,gene2,...,geneN --geneset-annot GenesetDatabase ）
6. 該位點(diǎn)最近的基因?qū)?yīng)的小鼠表型（--mouse-pheno）
7. 該位點(diǎn)最近的基因?qū)?yīng)的斑馬魚表型（--zebrafish-pheno）
8. 該位點(diǎn)最近的基因?qū)?yīng)的疾病表型（--ddd-annot）
9. 該位點(diǎn)最近的基因在pubmed數(shù)據(jù)庫中的相關(guān)文獻(xiàn)（--phenotype-term searchTerm and --pubmed-mining）

java -jar kggseq.jar --vcf-file path/to/file1 --ped-file path/to/file2 --phenotype-term search+Term1,search+Term2 --pubmed-mining

PubMedIDIdeogram : PubMed ID of articles in which the term and the cytogeneic position of the variant are co-mentioned
PubMedIDGene : PubMed ID of articles in which the term and the gene containing the variant are co-mentioned

10. 該位點(diǎn)對應(yīng)基因在OMIM數(shù)據(jù)庫中的疾病術(shù)語（ --omim-annot）
11. 該位點(diǎn)在COSMIC數(shù)據(jù)庫中的注釋信息（--cosmic-annot）

4. 位點(diǎn)功能的預(yù)測（Prediction at variant）

主要包括數(shù)據(jù)庫：i.e., SIFT, PolyPhen2, GWAVA and CADD

1. 預(yù)測孟德爾疾病致病位點(diǎn)（--db-score dbnsfp --mendel-causing-predict）

java -jar kggseq.jar --vcf-file path/to/file1 --ped-file path/to/file2 --db-score dbnsfp --mendel-causing-predict best --filter-nondisease-variant

該方法通過對位點(diǎn)進(jìn)行邏輯回歸建模，評估前10個AUC曲線來確定候選位點(diǎn)。
也可以指定數(shù)據(jù)庫進(jìn)行打分預(yù)測（--mendel-causing-predict 4,6,7）

2. 預(yù)測癌癥中的高頻體細(xì)胞突變（--db-score dbnsfp --cancer-mut-predict）
3. 預(yù)測復(fù)雜疾病的致病位點(diǎn)

• 對非編碼位點(diǎn)使用10個數(shù)據(jù)庫進(jìn)行打分預(yù)測（--db-score dbncfp_known [or dbncfp_all] --regulatory-causing-predict ）
  10個數(shù)據(jù)庫包括：SuRFR、GWAS3D、FunSeq2、FATHMM-MKL、FunSeq、DANN、CADD_CScore、 GWAVA_Unmatched、GWAVA_TSS、GWAVA_Region。
• 通過添加細(xì)胞類型特異性調(diào)控甲基化分子標(biāo)記預(yù)測致病基因（--db-score dbncfp_known [or dbncfp_all] --regulatory-causing-predict all --cell）
• 數(shù)據(jù)庫資源包括（1）16 ENCODE cell types，默認(rèn)采用GM12878細(xì)胞系。
  （2）127 RoadMap human reference epigenomes ，默認(rèn)采用E116 （GM12878 Lymphoblastoid Cell Line）

java -jar kggseq.jar --vcf-file path/to/file1 --db-score dbncfp_known --regulatory-causing-predict all --cell GM12878 
java -jar kggseq.jar --vcf-file path/to/file1 --db-score dbncfp_known --regulatory-causing-predict all --cell E116 

4. 基因特征特異性模型（--db-score dbncfp_known [or dbncfp_all] --ldgf-func-predict）
   用10個數(shù)據(jù)庫的打分預(yù)測非編碼區(qū)位點(diǎn)的致病情況

5. 針對基因的預(yù)測（Prediction at gene）

1. 基因的致病性預(yù)測（--patho-gene-predict）
   通過對Clinical Genomic Database數(shù)據(jù)庫的打分，得到基因的致病打分。

java -jar kggseq.jar --vcf-file path/to/file1 --ped-file path/to/file2 --out path/to/prefixname --excel --db-gene refgene --patho-gene-predict 

2. 表型挖掘（--phenotype-term searchTerm and --phenolyzer-prediction）
   通過使用phenolyzer工具對表型相關(guān)基因進(jìn)行排序，尋找到候選基因。

java -jar kggseq.jar --vcf-file path/to/file1 --ped-file path/to/file2 --out path/to/prefixname --excel --db-gene refgene --phenotype-term searchTerm1,searchTerm2 --phenolyzer-prediction

3. 組織特異性表達(dá)（--tissue-spec-annot）
   使用GTEx數(shù)據(jù)庫注釋基因的表達(dá)情況。

java -jar kggseq.jar --vcf-file path/to/file1 --ped-file path/to/file2 --out path/to/prefixname --excel --db-gene refgene --tissue-spec-annot 

6. 統(tǒng)計(jì)檢驗(yàn)

1. 癌癥驅(qū)動基因的負(fù)荷檢驗(yàn)-ITER或者WITER或者WITER-ref方法（--iter-gene-coding）

java -jar witer.jar --maf-file path/to/file1.maf --out path/to/prefixname [other tags] --iter-gene-coding
java -jar witer.jar --maf-file path/to/file1.maf --out path/to/prefixname [other tags] --witer-gene-coding
java -jar witer.jar --maf-file path/to/file1.maf --out path/to/prefixname [other tags] --witer-gene-coding  --ref-mut-file path/to/referenceDataset --cancer-label XXX

2. 復(fù)雜疾病易感基因罕見突變位點(diǎn)負(fù)荷檢驗(yàn)RUNER方法（--runer-gene-coding）

java -jar witer.jar --nt 6 [or more threads] --vcf-file path/to/file1 --ped-file path/to/file2 \
--out path/to/prefixname [or other tags] --excel --qqplot --gty-qual 20.0 --gty-sec-pl 20 --gty-dp 8 --gty-af-ref 0.05 \
--gty-af-het 0.25 --gty-af-alt 0.75 --ignore-indel --vcf-filter-in PASS --min-obs-rate 0.9 --hwe-all 0.001 \
--max-allele 3 --db-gene refgene,gencode  --db-filter gadexome.[an ancestry tag],gadgenome.[an ancestry tag],1kgeas201305 \
--rare-allele-freq 0.01   --gene-freq-score [an ancestry tag] \
--db-score dbnsfp --disease-causing-predict best --runer-gene-coding 

• 也可以選擇編碼區(qū)等位點(diǎn)進(jìn)行分析（ --response-gene-feature 2,3,4,5,6）
• 負(fù)荷檢驗(yàn)后可以挖掘最相關(guān)的文獻(xiàn)（--phenotype-term [disease name] --pubmed-mining-top-gene 10）
• 不將位點(diǎn)功能打分加入負(fù)荷檢驗(yàn)?zāi)Ｐ停?-uwruner-gene-coding）
• 使用本地control樣本的頻率進(jìn)行檢驗(yàn)（--gene-freq-score CONTROL）
• 加入其它協(xié)變量進(jìn)行檢驗(yàn)（ --gene-cov-file）,如：

Gene    score1  score1 ...
SAMD11  2.114   0.076914032 ...
NOC2L   2.345   0.046528595 ...
KLHL17  1.989   0.026329064 ...
PLEKHN1 2.2 0.047853652 ...
... ... ... ...

7. 對位點(diǎn)的關(guān)聯(lián)分析

1. 卡方檢驗(yàn)（--var-assoc）

java -jar ./kggseq.jar --buildver hg19  --var-assoc --p-value-cutoff 0.01 --multiple-testing benfdr --qqplot

2. 使用 Rvtests軟件的score, wald, exact, dominantExact, famLRT, famScore, famGrammarGamma 和 firth models等模型進(jìn)行檢驗(yàn)（ --rvtest-var）

java -jar ./kggseq.jar --rvtest-var score,wald,exact,... --rvtest-vcf --rvtest-remove-set \
--vcf-file path/to/vcf/file --ped-file path/to/pedigree/file --phe phenotypeName \
--cov covariateName1,covariateName2,...,covariateNameN \
--out path/to/prefixname --excel --db-gene refgene 

8. 對基因的關(guān)聯(lián)分析

1. 使用SKAT軟件以基因?yàn)閱挝贿M(jìn)行檢驗(yàn)（--skat-gene）
   包括SKAT,SKATO 和 Burden test算法。

java -jar ./kggseq.jar  --skat-gene --vcf-file path/to/vcf/file --ped-file path/to/pedigree/file \
--phe phenotypeName --cov covariateName1,covariateName2,...,covariateNameN \
--out path/to/prefixname --excel \
--db-gene refgene [--p-value-cutoff 0.01 --multiple-testing benfdr --qqplot]

2. 使用Rvtests（--rvtest-gene）軟件進(jìn)行檢驗(yàn)

java -jar ./kggseq.jar --rvtest-gene cmc,cmat,price,skat[nPerm=1000:alpha=0.001:beta1=1:beta2=20],... --rvtest-vcf \
--rvtest-remove-set --vcf-file path/to/vcf/file --ped-file path/to/pedigree/file --phe phenotypeName \
--cov covariateName1,covariateName2,...,covariateNameN \
--out path/to/prefixname --excel --db-gene refgene 

3. 使用SKAT軟件進(jìn)行基因通路分析（--skat-geneset）

java -jar ./kggseq.jar  --skat-geneset --geneset-file path/to/geneset/file --vcf-file path/to/vcf/file \
--ped-file path/to/pedigree/file --phe phenotypeName \
--cov covariateName1,covariateName2,...,covariateNameN \
--out path/to/prefixname --excel \
--db-gene refgene [--p-value-cutoff 0.01 --multiple-testing benfdr --qqplot]

基因通路需要存放在一個txt文件中，格式如下：
注：如果報(bào)錯，可能是因?yàn)槲募竺嬗锌招校堄?--lib-update命令更新KGGSeq到最新版，或直接刪除空行即可。

[Column 1: GeneSet ID]
[Column 2: GeneSet URL]
[Column 3..N: Gene Symbols in the geneset; separated by spaces].
No title row is required!!!
e.g.,
----------------------------------------------
KEGG_GLYCOLYSIS_GLUCONEOGENESIS http://www.broadinstitute.org/KEGG_GLYCOLYSIS_GLUCONEOGENESIS.html LDHC LDHB LDHA
KEGG_CITRATE_CYCLE_TCA_CYCLE http://www.broadinstitute.org/KEGG_CITRATE_CYCLE_TCA_CYCLE.html GJB1 OGDHL OGDH PDHB IDH3G LDHB
...
...

4. 使用RVTEST軟件進(jìn)行基因通路分析（--rvtest-geneset）
   最新版的 kggseq增加了參數(shù)--rvtest-min-var 2 可以去除位點(diǎn)小于2的基因。

java -jar ./kggseq.jar  --rvtest-geneset cmc,skato[nPerm=1000:alpha=0.001:beta1=1:beta2=20],... --rvtest-vcf \
--rvtest-remove-set --geneset-file path/to/geneset/file --vcf-file path/to/vcf/file \
--ped-file path/to/pedigree/file --phe phenotypeName --cov covariateName1,covariateName2,...,covariateNameN \
--out path/to/prefixname --excel \
--db-gene refgene [--p-value-cutoff 0.01 --multiple-testing benfdr --qqplot]

5. 基因富集分析（--geneset-enrichment-test）
   查看給定的基因是否在某些給定的基因通路上排在前列

java -jar kggseq.jar --vcf-file path/to/file1 --ped-file path/to/file2 --out path/to/prefixname [other tags] --geneset-enrichment-test
java -jar kggseq.jar --vcf-file path/to/file1 --ped-file path/to/file2 --out path/to/prefixname --db-filter 1kg201204,dbsnp137,ESP6500AA,ESP6500EA  --rare-allele-freq 0.01  --db-gene refgene --gene-mutation-rate-test --geneset-enrichment-test --geneset-db <GenesetDatabase>

可以自己設(shè)定通路基因（--geneset-file path/to/geneset/file）

9. 簡單的畫圖功能

1. MAF柱形圖（reference dataset --mafplot-ref or in sample --mafplot-sample）

java -jar kggseq.jar --vcf-file path/to/file1 --ped-file path/to/file2 --out path/to/prefixname --db-filter 1kg201204,dbsnp141,ESP6500AA,ESP6500EA  --allele-freq 0,1  --mafplot-ref
java -jar kggseq.jar --vcf-file path/to/file1 --ped-file path/to/file2 --out path/to/prefixname  --mafplot-sample

2. QQ圖

java -jar kggseq.jar --vcf-file path/to/file --out path/to/prefixname --qqplot --skat-gene [or] --var-assoc

10. 其他一些實(shí)用功能

1. 計(jì)算LD和基因型相關(guān)性（--calc-ld）

java -Xmx15g -jar kggseq.jar --vcf-file path/to/file --ped-file path/to/file --out path/to/prefixname \
--calc-ld --db-filter 1kg201204,dbsnp141,ESP6500AA,ESP6500EA  \
--allele-freq 0.05,0.95   --regions-in chr4:21212-233454

2. 刪除高LD的位點(diǎn)（--ld-prune X）

java -Xmx15g -jar kggseq.jar --vcf-file path/to/file \
--ped-file path/to/file --out path/to/prefixname \
--ld-prune 0.9 --filter-sample-maf-le 0.05 --o-svcf