肿瘤基因测序数据的常规分析大家应该是耳熟能详了，从上游测序fastq文件拿到每个病人的somatic突变信息，包括snv和cnv，无论是肿瘤全外显子和肿瘤基因panel测序的技术，都大同小异。我们生信入门马拉松授课团队的教学队长（nickier）在《生信技能树》和《生信菜鸟团》公众号都有相关知识整理专辑，比如《肿瘤外显子》专栏的目录（节选）如下：
- （一）读文献并且下载测序数据
- （二）质控与去接头
- （三）比对
- （四）比对结果的质控
- （番外篇）bam文件载入igv可视化
- （五）GATK的最佳实践
- （六）vcf文件的注释及ANNOVAR的使用
- （七）maftools可视化
- （八）不同注释软件的比较（上）：安装及使用
- （八）不同注释软件的比较（中）：注释后转成maf文件
- （八）不同注释软件的比较（下）：可视化比较maf文件
目前绝大部分医院的肿瘤相关突变数据的上游测序fastq文件处理都是公司代劳，流程是经典的GATK套装，主要是耗费计算机资源而已。
但是拿到了分析好的突变信息后续的理解又是另外一种难点了，比如针对肿瘤常见突变位点设计引物，我们随机选取3个位点来举例：
- KRAS基因的G13D突变
- AKT1基因的E17K突变
- EGFR基因的V843I突变
上面的肿瘤突变描述都是氨基酸突变，需要遵循人类基因组变异协会（HGVS：Human Genome Variation Society）（http://www.hgvs.org/mutnomen/）规则。主要是：
- 基因组参考序列（以前缀“g.”表示）
- 蛋白质参考序列（以前缀“p.”表示）
这些突变位点很容易在COSMIC数据库找到记录：
- https://cancer.sanger.ac.uk/cosmic/mutation/overview?id=99805418 （ GRCh38, 12:25245347,）
- https://cancer.sanger.ac.uk/cosmic/mutation/overview?id=158336406 （GRCh38, 14:104780214）
- https://cancer.sanger.ac.uk/cosmic/mutation/overview?id=99681646 （GRCh38, 7:55191776 ）
接下来就是问题转换，根据肿瘤突变位点，在对应的版本的参考基因组上面获取其上下游的碱基序列，足够长度可以去设计引物。
很简单的代码如下所示：
代码语言：javascript复制library(BSgenome.Hsapiens.UCSC.hg38)
# https://snystrom.github.io/memes-manual/reference/get_sequence.html
library(memes)
pos = c('chr12:25245347','chr14:104780214','chr7:55191776')
updown <- function(x){
chr=strsplit(x,":")[[1]][1]
p = as.numeric(strsplit(x,":")[[1]][2])
paste0(chr,":",p-20,'-',p 20)
}
lapply(pos, updown)
get_sequence(unlist(lapply(pos, updown)), BSgenome.Hsapiens.UCSC.hg38)
让我蛮意外的，这个memes获取到的突变位点上下游碱基序列居然还是彩色的。
突变位点上下游碱基序列
有了这些碱基信息，设计引物对绝大部分生命科学领域小伙伴来说就是小菜一碟啦。
我看了看这个get_sequence {memes} 函数，它是A light wrapper around Biostrings::getSeq to return named DNAStringSets, from input genomic coordinates.
但是没有看到其函数源代码，可能是需要去下载memes包的github源代码去解析了，看看它是如何加颜色的。
学徒作业
下载gwas catalog文件，比如我好久之前下载的 gwas_catalog_v1.0.2-associations_e105_r2021-12-21.tsv 文件，重要的信息，如下所示：
代码语言：javascript复制$ cat gwas_catalog_v1.0.2-associations_e105_r2021-12-21.tsv |awk -F"t" '{print $12,$13,$13 1,$21}'|head
CHR_ID CHR_POS 1 STRONGEST SNP-RISK ALLELE
9 82904977 82904978 rs7045089-A
20 49084487 49084488 rs2075677-A
3 43812828 43812829 rs6772840-A
5 152808169 152808170 rs4958581-T
3 158470022 158470023 rs6787172-T
20 32643466 32643467 rs6141769-C
10 3668695 3668696 rs10795061-T
1 20847936 20847937 rs17449243-T
5 104356621 104356622 rs10067176-A
针对你最感兴趣的疾病，拿到它在gwas catalog文件所显示的gwas位点的上下游碱基序列。
