too young, too simple

0%

全外显子测序数据处理

Posted on Views: Valine:
Symbols count in article: 8.7k Reading time ≈ 16 mins.

背景

老板的朋友的朋友被诊断为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的想法算是打消了。回到这个例子，退一步讲假设我们看到了这三百多个突变中有一个渐冻症相关，难道能对他误诊与否真起什么作用吗？

我还是得多学习一个。