背景
老板的朋友的朋友被诊断为ALS,也就是肌萎缩性脊髓侧索硬化症(也就是我们常说的渐冻人症),心有不甘,用二代测序测了自己的全外显子,想看一看自己有没有携带文献中已经报道的突变,万一是误诊呢?
咱们不谈这个神奇的想法,一个优秀的工具人只看这个事情该咋做。
Pipeline
这个东西基本的过程就是把测序得到的原始fastq比对到参考基因组上,然后找到和参考基因组不同的variant
, 使用网络上有的数据库注释这些variant 信息,看看跟文献报道的ALS相关突变(或者是任意)有无重合。目前很多公司都提供类似的数据分析服务,但是网上能搜到的结题报告效果很一般,就有一种我这活我也能做啊为啥找你的感觉。
大概的流程是一样的,我在本文中给出一个简单的实现,这显然不是最优的流程,但是区别不会太大:
基本上流程分为三个主要部分:
- 数据预处理(拿到bam文件)
- Call Variant(拿到vcf文件)
- Annotation (拿到眼睛看得懂的文件)
数据预处理
拿到bam文件这个过程每个实验室都有自己的脚本,我就在Deng-lab 的 chip-seq的流程上改了改:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
|
import os
def listdir_nohidden(path): for f in os.listdir(path): if not f.startswith('.'): yield f
SAMPLES = ["MiaoLL"] fastaIndex="/home/skelviper/Desktop/ref/GRCh37d5/hs37d5.fa.fai" bowtieIndex="/home/skelviper/Desktop/ref/GRCh37d5/bt2Index/hs37d5"
rule all: input: expand("processed/bw/{sample}.bw",sample=SAMPLES), expand("processed/qc/{sample}qc",sample=SAMPLES), expand("processed/bw/{sample}.bw",sample=SAMPLES), shell: """ set +u; source activate; conda activate chip; set -u
multiqc ./processed set +u;conda deactivate; set -u """ rule fastqc: input: fq1="Rawdata/{sample}.R1.fq.gz", fq2="Rawdata/{sample}.R2.fq.gz", output: qcRes = directory("processed/qc/{sample}qc") threads: 10 shell: """ set +u; source activate; conda activate chip; set -u
mkdir -p processed/qc/{wildcards.sample}qc fastqc {input.fq1} {input.fq2} -o {output} -t {threads} set +u;conda deactivate; set -u """ rule fq2trim: input: fq1=rules.fastqc.input.fq1, fq2=rules.fastqc.input.fq2, output: trim1 = temp("processed/trim/{sample}/{sample}.R1_val_1.fq.gz"), trim2 = temp("processed/trim/{sample}/{sample}.R2_val_2.fq.gz"), trimDir = directory("processed/trim/{sample}/"),
threads: 10 shell: """ set +u; source activate; conda activate trim; set -u
trim_galore --cores {threads} --paired {input.fq1} {input.fq2} -o {output.trimDir} --fastqc set +u;conda deactivate; set -u """
rule mapping: input: trim1 = rules.fq2trim.output.trim1, trim2 = rules.fq2trim.output.trim2, output: sam = "processed/sam/{sample}.sam", threads: 10 log: logfile="processed/mappingStat/{sample}.sam.stat", shell:""" set +u; source activate; conda activate chip; set -u
mkdir -p processed/mappingStat bowtie2 -t -q -p {threads} -N 1 -L 25 -X 700 --no-unal -x {bowtieIndex} -1 {input.trim1} -2 {input.trim2} -S {output.sam} 2>> {log.logfile} set +u;conda deactivate; set -u """
rule generateDedupBam: input: sam = rules.mapping.output.sam, output: sortedBam = temp("processed/bamFile/{sample}.sorted.bam"), rmdupBam = "processed/bamFile/{sample}.sorted.rmdup.bam", marked_dup_metrics = "processed/bamFile/{sample}.marked_dup_metrics.txt"
shell:""" set +u; source activate; conda activate chip; set -u
samtools view -bt {fastaIndex} {input.sam} | samtools sort > {output.sortedBam} picard MarkDuplicates -I {output.sortedBam} -O {output.rmdupBam} -M {output.marked_dup_metrics} --REMOVE_DUPLICATES true --ASSUME_SORTED true samtools index {output.rmdupBam} set +u;conda deactivate; set -u """
rule readsVis: input: rmdupBam = rules.generateDedupBam.output.rmdupBam, output: bigwig = "processed/bw/{sample}.bw" shell:""" set +u; source activate; conda activate macs2; set -u
bamCoverage --normalizeUsing RPKM -b {input.rmdupBam} -o {output.bigwig} set +u;conda deactivate; set -u """
|
(回头应该在这里做个代码折叠啥的,有一点点影响阅读体验。
我的觉得最后的bam文件在任何的流程中都应该保留下来,这样你对下游不自信的时候,还总是可以返回去在igv里手动看一看bam上的比对情况。特别是对于这种Call variant 的项目,谁知道它call的对不对呢,随便看几个心里能有个底。
Call Variant
一般这个部分用的是gatk
, 一般大家说起call variant都会想起它。但是我这里用了一个叫DeepVariant
的软件, 根据一些文献报道,DeepVariant 在准确性和性能上都要比gatk强一点。
https://www.nature.com/articles/s41598-020-77218-4 (Fig4,DeepVariant流程的时间大幅小于GATK流程)
https://www.nature.com/articles/s41598-018-36177-7 (摘要,使用F-score 度量,DeepVariant在对indel 的处理上强于gatk,在snv的处理上基本一致)
DeepVariant 非常好用的一点在于官方提供了Docker环境,另外听说gatk有一些多线程不够多的问题(这一点我没有求证)。
1
| sudo docker run -v "${work_dir}":"/input" -v "${work_dir}":"/output" google/deepvariant:"1.1.0" /opt/deepvariant/bin/run_deepvariant --model_type=WES --ref=/input/hs37d5.fa --reads=/input/MiaoLL.sorted.rmdup.bam --output_vcf=/output/MiaoLL.vcf --output_gvcf=/output/MiaoLL.gvcf --num_shards=15 --logging_dir=/output/logs
|
相信的文档可以在deepvariant 的github页面上找到:
https://github.com/google/deepvariant
需要注意的是,docker虽然给你省了点配环境的时间,但是容器内的应用访问不了容器外的文件,这意味着如果你把iwork_dir里面的fa和bam都放的是链接,deepvariant会报错。不过写个snakemake copy一份文件过去,最后把结果mv出来,然后整个workdir当作temp也不是什么困难的事情。这样我们就得到了vcf文件。
VCF annotation
首先,上一步的vcf需要一些过滤,比如我将QUAL < 20 的行都去掉了。
1
| cat MiaoLL.vcf | awk '$6>20 {print $0}' > MiaoLL.Q20.vcf
|
vcf 的注释有很多种选择,一种比较简单的是使用vep。我懒得配环境,直接使用了在线版本,这个步骤需要注意的一点是要选择和上游mapping相同的参考基因组(38还是37)
我使用的等效命令行参数如下:
1
| ./vep --af --appris --biotype --buffer_size 500 --check_existing --distance 5000 --filter_common --mane --polyphen b --pubmed --regulatory --sift b --species homo_sapiens --symbol --transcript_version --tsl --cache --input_file [input_data] --output_file [output_file] --port 3337
|
拿到了这样一个结果:

把这个数据下载下来,你就有了很全面的注释信息。注意这里拿到了很多内含子突变,这是因为WES对DNA的富集不是精确的,总是有延伸出来的到内含子区域的部分。
Annotation
vep得到的原始注释信息还是太多了,并且其中很多是所谓的Benign,这些是我不关心的。稍微过滤一下不属于Benign 分类的位点:
可以看一下这里面有啥类型:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
| - affects affects,association association association,benign association,benign/likely_benign,pathogenic,drug_response,conflicting_interpretations_of_pathogenicity association,likely_benign association,not_provided,drug_response benign benign,affects benign,association benign,benign/likely_benign benign,benign/likely_benign,benign,_other benign,benign/likely_benign,likely_benign benign,benign/likely_benign,not_provided benign,conflicting_interpretations_of_pathogenicity benign,drug_response benign,drug_response,likely_benign benign,likely_benign benign/likely_benign benign/likely_benign,benign ……………省略一大堆
|
选一下我们可能会关心的:
1
| cat MiaoLL.Q20.vepRES.txt | awk '$57 != "-"{print $0}' | grep -e 'risk' -e 'path' -e 'drug' | grep -P 'rs[0-9]*' -o | sort | uniq > rsnums.txt
|
rs号对应的疾病可以在ncbi 的dbSNP上查到,但是dbSNP没有提供一个批量查询的接口。
好在这个网站非常简单,地址栏直接替换一下rs号就能搜索,结果就是一个表格。

我们可以做个小爬虫:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
|
import requests,random import csv from bs4 import BeautifulSoup import pandas as pd import time import argparse
parser = argparse.ArgumentParser(description="Crawler for rs-num to Disease name") parser.add_argument( "-i","--input", dest="rsnumsPath", type=str, action="store", help="input folder path of contactMatrix" ) parser.add_argument( "-o","--output", dest="outputResPath", type=str, action="store", help="input folder path of contactMatrix" ) args = parser.parse_args()
rsnums = [] with open(args.rsnumsPath) as f: rsnums = f.read().splitlines() def get_one_rs_table(rsnum): user_agents=['Mozilla/5.0 (Windows NT 6.1; rv:2.0.1) Gecko/20100101 Firefox/4.0.1', 'Mozilla/5.0 (Windows; U; Windows NT 6.1; en-us) AppleWebKit/534.50 (KHTML, like Gecko) Version/5.1 Safari/534.50', 'Opera/9.80 (Windows NT 6.1; U; en) Presto/2.8.131 Version/11.11'] headers={'User-Agent':random.choice(user_agents)} req=requests.get("".join(["https://www.ncbi.nlm.nih.gov/snp/",rsnum,"#clinical_significance"]),headers=headers) html=req.text bf=BeautifulSoup(html,features="lxml") bf_raw = bf.select("#clinical_significance_datatable")[0].prettify() df = pd.read_html(bf_raw)[0] df.insert(0,'rsNum',rsnum) time.sleep(1) return df
res = pd.DataFrame() for rsnum in rsnums: try: res = res.append(get_one_rs_table(rsnum)) except Exception as e: pass continue
res.columns = ['rsNum', 'ClinVar_Accession',"Disease_Names","Clinical_Significance"] res = res.query('Clinical_Significance != "Benign"') res = res.query('Clinical_Significance != "Not-Provided"') res = res.reset_index(drop="index")
res.to_csv(args.outputResPath,sep="\t")
|
运行一下:
1
| python ~/myData/project/202103/als_snp/scripts/rsnum2Diseases.py -i rsnums.txt -o tableres.csv
|
最后能得到一个类似这样的结果

三百多行,大概就是一个人眼能看的长度了。如果你想更进一步,大概可以做个词云或者NLP啥的。我不想把这个工作延伸到那么远。
讨论
在这里我们并没有看到ALS相关的突变,也许有些位点新发现没有收录到数据库中。不过还是找到了不少其他的。仔细看一下这些Risk Factor,这位朋友从高血压心脏病糖尿病到阿尔茨海默症一应俱全。我不排除这个人是在运气太差了的可能性,但更大可能性的,这些位点的存在本身什么都说明不了,同时易感所有疾病就是对所有疾病都不易感。我很好奇公司会怎么处理这样的报告,但是目前来看我自己测个WES的想法算是打消了。回到这个例子,退一步讲假设我们看到了这三百多个突变中有一个渐冻症相关,难道能对他误诊与否真起什么作用吗?
我还是得多学习一个。