kggseq软件用于测序后下游分析

1. 质控(Quality Control)

  1. 对基因型的质控(Genotype QC)
    参数如下:


  2. 对位点的质控(Variant QC)
    参数如下:


    image.png

    image.png
  3. 对样本的质控(Sample QC)
    可以使用kggseq软件把vcf转换为plink格式进行质控,包括 (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参数,后面加上数字,数字对应意义如下:


    image.png
  • 使用--double-hit-gene-trio-filter 或 --double-hit-gene-phased-filter参数提取复合杂合或者隐形突变位点
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)
  1. 使用基因特征过滤(Gene feature filtering)
    参数:--db-gene & --gene-feature-in (两个参数分别选择数据库和基因特征)
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
  1. 使用公共数据库的基因型频率进行过滤(--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

也可以对本地数据进行频率过滤(--local-filter-hard)

  1. 过滤掉重复区域(--superdup-filter)或者注释这些区域(--superdup-annot)
java -jar ./kggseq.jar --buildver hg19 --vcf-file path/to/file1 --ped-file path/to/file2 --superdup-filter
  1. 过滤掉indel,某些基因或某些区域
  • --ignore-indel
  • --ignore-snv
  • --genes-out
  • --regions-out
  1. 过滤掉假基因等其他因素造成的变异位点过多的基因(--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
  1. 过滤掉实验组/对照组独有的位点(--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. 选择性剪切位点(--scsnv-annot)
  3. 该位点是否位于CNV区域(--dgv-cnv-annot)
  4. 该位点是否与感兴趣的基因在同一蛋白通路上(--candi-list gene1,gene2,...,geneN --ppi-annot ppiDatabase )
  5. 该位点是否与感兴趣的基因在同一基因通路上(--candi-list gene1,gene2,...,geneN --geneset-annot GenesetDatabase )
  6. 该位点最近的基因对应的小鼠表型(--mouse-pheno)
  7. 该位点最近的基因对应的斑马鱼表型(--zebrafish-pheno)
  8. 该位点最近的基因对应的疾病表型(--ddd-annot)
  9. 该位点最近的基因在pubmed数据库中的相关文献(--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
  1. 该位点对应基因在OMIM数据库中的疾病术语( --omim-annot)
  2. 该位点在COSMIC数据库中的注释信息(--cosmic-annot)

4. 位点功能的预测(Prediction at variant)

主要包括数据库:i.e., SIFT, PolyPhen2, GWAVA and CADD

  1. 预测孟德尔疾病致病位点(--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

该方法通过对位点进行逻辑回归建模,评估前10个AUC曲线来确定候选位点。
也可以指定数据库进行打分预测(--mendel-causing-predict 4,6,7)

  1. 预测癌症中的高频体细胞突变(--db-score dbnsfp --cancer-mut-predict)
  2. 预测复杂疾病的致病位点
  • 对非编码位点使用10个数据库进行打分预测(--db-score dbncfp_known [or dbncfp_all] --regulatory-causing-predict )
    10个数据库包括:SuRFR、GWAS3D、FunSeq2、FATHMM-MKL、FunSeq、DANN、CADD_CScore、 GWAVA_Unmatched、GWAVA_TSS、GWAVA_Region。
  • 通过添加细胞类型特异性调控甲基化分子标记预测致病基因(--db-score dbncfp_known [or dbncfp_all] --regulatory-causing-predict all --cell)
  • 数据库资源包括(1)16 ENCODE cell types,默认采用GM12878细胞系。
    (2)127 RoadMap human reference epigenomes ,默认采用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 
  1. 基因特征特异性模型(--db-score dbncfp_known [or dbncfp_all] --ldgf-func-predict)
    用10个数据库的打分预测非编码区位点的致病情况

5. 针对基因的预测(Prediction at gene)

  1. 基因的致病性预测(--patho-gene-predict)
    通过对Clinical Genomic Database数据库的打分,得到基因的致病打分。
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 
  1. 表型挖掘(--phenotype-term searchTerm and --phenolyzer-prediction)
    通过使用phenolyzer工具对表型相关基因进行排序,寻找到候选基因。
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
  1. 组织特异性表达(--tissue-spec-annot)
    使用GTEx数据库注释基因的表达情况。
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. 统计检验

  1. 癌症驱动基因的负荷检验-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
  1. 复杂疾病易感基因罕见突变位点负荷检验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 
  • 也可以选择编码区等位点进行分析( --response-gene-feature 2,3,4,5,6)
  • 负荷检验后可以挖掘最相关的文献(--phenotype-term [disease name] --pubmed-mining-top-gene 10)
  • 不将位点功能打分加入负荷检验模型(--uwruner-gene-coding)
  • 使用本地control样本的频率进行检验(--gene-freq-score CONTROL)
  • 加入其它协变量进行检验( --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. 对位点的关联分析

  1. 卡方检验(--var-assoc)
java -jar ./kggseq.jar --buildver hg19  --var-assoc --p-value-cutoff 0.01 --multiple-testing benfdr --qqplot
  1. 使用 Rvtests软件的score, wald, exact, dominantExact, famLRT, famScore, famGrammarGamma 和 firth models等模型进行检验( --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. 对基因的关联分析

  1. 使用SKAT软件以基因为单位进行检验(--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]
  1. 使用Rvtests(--rvtest-gene)软件进行检验
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 
  1. 使用SKAT软件进行基因通路分析(--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文件中,格式如下:
注:如果报错,可能是因为文件后面有空行,请用 --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
...
...
  1. 使用RVTEST软件进行基因通路分析(--rvtest-geneset)
    最新版的 kggseq增加了参数--rvtest-min-var 2 可以去除位点小于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]

  1. 基因富集分析(--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>

可以自己设定通路基因(--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
  1. QQ图
java -jar kggseq.jar --vcf-file path/to/file --out path/to/prefixname --qqplot --skat-gene [or] --var-assoc

10. 其他一些实用功能

  1. 计算LD和基因型相关性(--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
  1. 删除高LD的位点(--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 
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 204,684评论 6 478
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 87,143评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,214评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,788评论 1 277
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,796评论 5 368
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,665评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,027评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,679评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 41,346评论 1 299
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,664评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,766评论 1 331
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,412评论 4 321
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,015评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,974评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,203评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,073评论 2 350
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,501评论 2 343

推荐阅读更多精彩内容

  • 刘小泽写于18.9.22今天先更新一半🤤不管是转录组,还是芯片数据,或者其他有关基因的组学分析,每当数据分析到后面...
    刘小泽阅读 5,732评论 4 25
  • 一.对今日学人力资源管理的总结: 打招聘电话及面试技巧: 1.接通电话后问是否方便接电话,工作是否落实 2.邀请来...
    上海细胞治疗小雅阅读 155评论 0 0
  • 有多少物,不舍得丢弃,暂时存放于某处,却在岁月的长河中渐渐被遗忘。被遗忘还是好的,若哪一天又被翻到,仍是带来苦恼,...
    甘棠夜话阅读 232评论 0 0
  • 为什么会有桥? 路走到头了,没有路便有了桥。 听说,人死后也要走过一座奈何桥,忘记尘世的所有记忆,那一辈子翻篇了,...
    夏若禅阅读 535评论 0 3
  • 人这一生啊,总要坚持一些东西 坚持做一些事情 坚持喜欢一些人 尽管有时候 这些坚持看起来有点傻 可是谁没有傻傻的时...
    纳西奇麒阅读 211评论 0 0