0%

全外显子测序数据处理

背景

老板的朋友的朋友被诊断为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
##########################################
# Standard ChIPseq/ATAC pipeline #
# @author zliu #
# @date 2021/03/04 #
##########################################

import os

def listdir_nohidden(path):
for f in os.listdir(path):
if not f.startswith('.'):
yield f


############CONFIG###########
#SAMPLES = list(listdir_nohidden("./Rawdata"))
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),
#expand("processed/narrowpeaks/{sample}.narrowpeak",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

拿到了这样一个结果:

image-20210314121005011

把这个数据下载下来,你就有了很全面的注释信息。注意这里拿到了很多内含子突变,这是因为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号就能搜索,结果就是一个表格。

image-20210314183831865

我们可以做个小爬虫:

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
# Crawler for converting rs-num to Disease name, data from dbSNP,NCBI.
# @author zliu 2021/3/13
# Usage: python rsnum2Diseases.py -i input.rsnum.txt -o output.res.tsv


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

最后能得到一个类似这样的结果

Result

三百多行,大概就是一个人眼能看的长度了。如果你想更进一步,大概可以做个词云或者NLP啥的。我不想把这个工作延伸到那么远。

讨论

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

我还是得多学习一个。