Predicting effects of noncoding variants with deep learning–based sequence model | 基於深度學習的序列模型預測非編碼區

Predicting effects of noncoding variants with deep learning–based sequence model 

PDF

Interpreting noncoding variants - 很是好的學習資料

這篇文章的第一個亮點就是直接從序列開始分析，第二就是使用深度學習得到了很好的預測效果。 

This is, to our knowledge, the first approach for prioritization of functional variants using de novo regulatory sequence information.

該寫的文章裏都寫得很是清楚了，再看不懂就是本身的問題了。

 

 

理解數據

List of all publicly available chromatin feature profile files used for training DeepSEA

ENCODE有哪些public的數據？本文中用了哪些數據來訓練？

• DHS - wgEncodeAwgDnaseUniform
• TF binding - wgEncodeAwgTfbsUniform
• Histone marks - roadmap

To prepare the input for the deep convolutional network model, we split the genome into 200-bp bins. For each bin we computed the label for all 919 chromatin features; a chromatin feature was labeled 1 if more than half of the 200-bp bin is in the peak region and 0 otherwise.

這919 chromatin features都是些什麼，怎麼來的？(125 DNase features, 690 TF features, 104 histone features)

The 1,000-bp DNA sequence is represented by a 1,000 × 4 binary matrix, with columns corresponding to A, G, C and T. The 400-bp flanking regions at the two sides provide extra contextual information to the model.

訓練集的上下游是加了400bp，可是訓練的時候如何利用這額外的序列呢？

這裏的input就很明確了：1000x4xN，而後每個1000x4都對應了919個結果，也就是input的維度爲1000，predict出來結果的維度爲919.

這樣，給定任何一個長度爲1000的序列，就能獲得一個919維度的關於chromatin feature的結果。

關於第一個model的實現就基本完成了。

 

那麼第二個model呢？

For positive standards we used single-nucleotide substitution variants annotated as regulatory mutations in the HGMD professional version 2014.4 (ref. 17), 

eQTL data from the GRASP 2.0.0.0 database with a P-value cutoff of 1 × 10−10 (ref. 1)

GWAS SNPs downloaded from the NHGRI GWAS Catalog on 17 July 2014 (ref. 18).

額外：Coding variants were filtered on the basis of the UCSC build hg19 knownGene track22.

SNP和eQTL的數據是如何被結合到當前的model裏的？

 

To compute predicted chromatin effects of variants using the DeepSEA model, for each SNP, we obtained the 1,000-bp sequence centered on that variant based on the reference genome (specifically, the sequence is chosen so that the variant was located at the 500th nucleotide). Then we constructed a pair of sequences carrying either the reference or alternative allele at the variant position.

目的就是要預測每個variant對919個chromatin feature的影響，此時咱們就知道了任意一個序列的其中任意一個位點的SNP致使的919個chromatin feature差別結果。

 

而後這個如何與eQTL結合呢？

For each of the three variant types, HGMD single-nucleotide substitution regulatory mutations, eQTLs and GWAS SNPs, we trained a regularized logistic regression model, using the XGBoost implementation

這個分類模型是用來預測什麼的呢？

The HGMD regulatory mutation model was trained with L1 regularization parameter 20 and L2 regularization parameter 2,000 for ten iterations. eQTL and GWAS SNP models were trained with L1 regularization parameter 0 and L2 regularization parameter 10 for 100 iterations.

能夠看到這兩個數據是分開train的，這些數據裏面分別都有些什麼信息？這個model到底想用什麼來訓練什麼？SNP來訓練出eQTL嗎？

HGMD：Substitutions causing regulatory abnormalities are logged in with thirty nucleotides flanking the site of the mutation on both sides. The location of the mutation relative to the transcriptional initiation site, initiaton codon, polyadenylation site or termination codon is given.

GRASP: Genome-Wide Repository of Associations Between SNPs and Phenotypes 這裏的eQTL應該是每個SNP高度相關的基因表達的Rank。

 

第二個model介紹比較簡單，由於不少人都在作這方面的工做，文章裏也列出了不少個比較的工具。

但由於我沒作過，因此暫時很難理解這部分的工做。

 

理解問題

無問題，不學習：

只關注non-coding區，那coding區的就無論了嗎？coding區是直接對蛋白產生影響，同時還會對某些lncRNA的產生影響。須要紮實的基因組學知識。

預測variants的effect，這裏的effect怎麼直觀理解？基因表達是effect嗎？

這是一個什麼樣的深度學習模型？input的feature是什麼？output的result是什麼？

用了什麼數據來訓練這個模型？ENCODE和Roadmap

directly learns a regulatory sequence code from large-scale chromatin-profiling data, enabling prediction of chromatin effects of sequence alterations with single-nucleotide sensitivity. 

該模型如何考慮了eQTL和疾病相關的variants？

關於non-coding variant你有什麼理解？數量大，功能不明確，經過eQTL咱們確實發現了一些與基因表達密切相關的位點。但顯然，某個基因會有不少個eQTL位點，這些位點是如何具體影響到某個基因的表達的，咱們不得而知。如今已經發現了eQTL位點在不一樣組織裏具備不一樣的效應，說明了基因的選擇性表達。因爲中心法則的限制，變異能影響的有轉錄、轉錄調控和翻譯。目前關於基因組序列的運做功能咱們知之甚少，但它又是一切的源泉，因此最基本的變異帶來的影響其實很複雜的，30億個鹼基，可能的變化組合太多了。從細節開始分析這些variant的功能的策略顯然是行不通的。

 

咱們的理解和認知極限，爲何咱們對複雜系統的分析是如此的艱難？一個真實存在並高效運轉的系統，咱們卻猜不透摸不着。

一個變異很好理解，找出細胞內全部與該位點直接相關的生物過程，再偵測後續如何發生連鎖反應，最終會致使寫什麼結果！

基因組的變異是時刻在發生的，體細胞、生殖細胞，咱們真正面對的是成千上萬個突變，這些突變綜合在一塊兒是如何相互影響，致使最終的表型的變化的？這複雜度可就呈指數上升了，中間的過程能夠很複雜很瑣碎，咱們人真正關注的仍是最終的表型，生死和慾望。

 

no method has been demonstrated to predict with single-nucleotide sensitivity the effects of noncoding variants on transcription factor (TF) binding, DNA accessibility and histone marks of sequences.

A quantitative model accurately estimating binding of chromatin proteins and histone marks from DNA sequence with singlenucleotide sensitivity is key to this challenge.

 

Therefore, accurate sequence-based prediction of chromatin features requires a flexible quantitative model capable of modeling such complex dependencies—and those predictions may then be used to estimate functional effects of noncoding variants.基於序列的，能夠考慮到更多複雜的狀況，這比單純的依賴現有的motif要更加準確。

We first directly learn regulatory sequence code from genomic sequence by learning to simultaneously predict large-scale chromatin-profiling data, including TF binding, DNase I sensitivity and histone-mark profiles. 

 

sequence context

multilayer

share learned predictive sequence features

 

 

通俗來講：

這篇文章裏有兩個model，第一個是根據序列的特徵來預測enhancer和promoter，而後再基於enhancer和promoter的epigenome的內容來預測eQTL。

model1：TF binding, DNase I sensitivity and histone-mark profiles

 

文章晦澀難懂能夠直接先跑下代碼試下：

代碼在這：http://deepsea.princeton.edu/job/analysis/create/

三種可選的input文件：

FASTA

>[known_CEBP_binding_increase]GtoT__chr1_109817091_109818090
GTGCCTCTGGGAGGAGAGGGACTCCTGGGGGGCCTGCCCCTCATACGCCATCACCAAAAGGAAAGGACAAAGCCACACGC
AGCCAGGGCTTCACACCCTTCAGGCTGCACCCGGGCAGGCCTCAGAACGGTGAGGGGCCAGGGCAAAGGGTGTGCCTCGT
CCTGCCCGCACTGCCTCTCCCAGGAACTGGAAAAGCCCTGTCCGGTGAGGGGGCAGAAGGACTCAGCGCCCCTGGACCCC
CAAATGCTGCATGAACACATTTTCAGGGGAGCCTGTGCCCCCAGGCGGGGGTCGGGCAGCCCCAGCCCCTCTCCTTTTCC
TGGACTCTGGCCGTGCGCGGCAGCCCAGGTGTTTGCTCAGTTGCTGACCCAAAAGTGCTTCATTTTTCGTGCCCGCCCCG
CGCCCCGGGCAGGCCAGTCATGTGTTAAGTTGCGCTTCTTTGCTGTGATGTGGGTGGGGGAGGAAGAGTAAACACAGTGC
TGGCTCGGCTGCCCTGAGGTTGCTCAATCAAGCACAGGTTTCAAGTCTGGGTTCTGGTGTCCACTCACCCACCCCACCCC
CCAAAATCAGACAAATGCTACTTTGTCTAACCTGCTGTGGCCTCTGAGACATGTTCTATTTTTAACCCCTTCTTGGAATT
GGCTCTCTTCTTCAAAGGACCAGGTCCTGTTCCTCTTTCTCCCCGACTCCACCCCAGCTCCCTGTGAAGAGAGAGTTAAT
ATATTTGTTTTATTTATTTGCTTTTTGTGTTGGGATGGGTTCGTGTCCAGTCCCGGGGGTCTGATATGGCCATCACAGGC
TGGGTGTTCCCAGCAGCCCTGGCTTGGGGGCTTGACGCCCTTCCCCTTGCCCCAGGCCATCATCTCCCCACCTCTCCTCC
CCTCTCCTCAGTTTTGCCGACTGCTTTTCATCTGAGTCACCATTTACTCCAAGCATGTATTCCAGACTTGTCACTGACTT
TCCTTCTGGAGCAGGTGGCTAGAAAAAGAGGCTGTGGGCA
>[known_FOXA2_binding_decrease]ReferenceAllele__chr10_23507864_23508863
CTTCTTTTTATCTCTTAACTAACTTACAATTTCTTACGTGATTTTAAAACTTGTTTTTCTATTTAAAACAACAGGGGCAA
CTGAACTTCACTTTCAAACAATATTTATTTCTATAAATCAGTGCAAAACATACTTATTGAAAATATATCTTGGGTCCAAG
GCTTCAAAGGGTAAAAAGAAAGATTTTAAATTATATCTAATATGTTACAATTGTTCTGTCCTTTAAAAACCTTTTCAGAT
CACCCCCTGGATGATTCTTCCCTAGAAGTCTCAGAGAATTAACAACACAATGTAATCTAGGTTTAAATTTGGGTTTCTCC
TGTGTTTCAGATACTGATGTTTGAGCTTTCTCTTCCTGACAAGCCACTTAAAGAGTCACTGTTACTTTGAGGTTTTATCT
GTAAGATTCGTGTCTTTTGGGCTCATTAAGAACATTTCCAAAGATTACAATGTCAATAGCACCTAATTACTGGACTGTGA
GAAAGGTCTTCTTGAGTACATAAAATCTGTGGCAGTGCACAGTACACAATGGGCAGCTCAGATCCCAAATTTTATCACAA
GTAAGTAGCAAACAAATTAATAATGTTACCTGTGCTCTCTTGGATAATTACTACTGCATAAAAACTGCTTTGAAATGTTG
CAGATAGTATTGTACCTCATTTTTTTAATCCCCTTAGAGTAACAAGGATTTATTTGTCTCAAACTTTCTATGTTGCATGC
ACCACTTGACTTTCTTGTTCTGTTTAGAATTTTTAGAACTTGCAACATAACAAAAAATCATTTTTAACCAGCCTAGGAAG
GACATATCACCTGATGTAACATTATTTTAAATTATATTTTGTATTTTACTTTACTCTTTTCAAAACATATACTGTATGTT
TTGATACTATTGCTAGATTTTATTTTTTACTTATGCCTGGTAGAAAATCAGCTATTAAAGAAGCAGAGGAGGCTGGACAC
AGTGGTTCATGTCTGTAATCGCTAGCACTTTGAAAGAGTA
>[known_GATA1_binding_increase]TtoC__chr16_209210_210209
GGGCTTAGACAGAGGAGGGGAGGATTCAGATTTTAAATGGGTTGGCCACTGTAGGTCTATTAACGTGGTGACATTTGAGG
GAGTGGCAATACTAGGGAAGGGGCTTCAGGGGAGTGGCCAGGAGCTAGGGATAGAGGGAGGGAGGACAGGAGGCCTTGTC
TGTCTTTTCCTCCATATGTAAGTTTCAGGAGTGAGTGGGGGGTGTCGAGGGTGCTGTGCTCTCCGGCCTGAGCCTCAGGA
AGGAAGGGCAGTAGTCAGGGATGCCAGGGAAGGACAGTGGAGTAGGCTTTGTGGGGAACTTCACGGTTCCATTGTTGAGA
TGATTTGCTGGAGACACACAGATGAGGACATCAAATACATCCCTGGATCAGGCCCTGGGGCCTGAGTCCGGAAGAGAGGT
CTGTATGGACACACCCATCAATGGGAGCACCAGGACACAGATGGAGGCTAATGTCATGTTGTAGACAGGATGGGTGCTGA
GCTGCCACACCCACATTATCAGAAAATAACAGCACAGGCTTGGGGTGGAGGCGGGACACAAGACTAGCCAGAAGGAGAAA
GAAAGGTGAAAAGCTGTTGGTGCAAGGAAGCTCTTGGTATTTCCAATGGCTTGGGCACAGGCTGTGAGGGTGCCTGGGAC
GGCTTGTGGGGCACAGGCTGCAAGAGGTGCCCAGGACGGCTTGTGGGGCACAGGTTGTGAGAGGTGCCCTGGACGGCTTG
TGGGGCACAGGCTGTGAGAGGTGCCCAGGACGGCTTGTGGGGCACAGGCTGTGAGGGTGCCCGGGACGGCTTGTGGGGCA
CAGGTTGTGAGAGGTGCCCGGGACGGCTTGTGGGGCACAGGTTTCAGAGGTGCCCGGGACGGCTTGTGGGGCACAGGTTG
TGAGAGGTGCCCGGGACGGCTTGTGGGACACAGGTTGTGAGAGGTGCCTGGGACGGCTTGTGGGGCACAGGCTGTGAGGG
TGCCTGGGACGGCTTGTGGGGCACAGGTTGTGAGAGGTGC
>[known_FOXA1_binding_increase]CtoT_chr16_52598689_52599688
GGCTCAAGCAGTCCTCCCATCTAGGCTTCCCAAAATGCTGGGATTACAGACATGAGCCACTGCACCCAGCCACAAAGATA
ACCTAAAGATGTGTTTACTTTGACCCAGGCAGTAGTTTAAAAAAGTTTTAATTTGTTGTTCACATTTAAAAACTGGACAA
TTTCTACATAAAAATCTGAATTACTCATGTCTCTTAAAAAAATAACATCTAGCAATGGTAGGCCCACATTCCTTCCTGAA
AATAATTAGCTGGGAAAGAGTAGGGACTGACCCCTTTAGACACGGTATAAATAGCATGGGAGTTGATCAGTAAATATTTG
CTGAATGAAAGAATACATGAATGAAAAGTCAGAGCCCTATAGGTCAGCATGGACGGCGGTAAAGGAACCTGGCTGAGCCT
GAAAGAGAATGTGATCTAAGATTAAATCCAGGATATGCTGGTAAATGTTTAACAGCCAACTCTTTGGGGAGGAAAAAAGT
CCCAATTTGTAGTGTTTGCTGATTATTGTGATGTAAATACTCCCATCATGACCAATTTCAAGCTACCAACATGCTGACAC
TGAACTTGGAGTTGGAAGGAGATGAACAGGCATAATCAGGTCTCGTGAGATGGCCCAAGCCGGCCCCAGCACTCCACTGT
TATATATGAGGCTAGAATTACTACATAACTGGAATAGCAACTTTCTGGACCATATGCCTGGAACACAGCAGGTGCTGAAT
AAATGTTTGTTGATCCAGGAACTGACTGTGTTGAAGCCCACAGATGGGAAATCAGTAGAAGGCAGGTAAGAGTAAAAAGA
AGGGCAGAGAATTGGGGGTACAGACCCCTGAACCATAAGTCAGAGGAATGTTGTACATGTTTTCAGATCCCTCACTGGTC
AAATGAAGGCAAAGGGTTAGATCTCTCCAAATCTTTAGAGGGACATGATGTAACTCCATTAAGTAACTCAGTGATTTTCA
ACATTAAAAAGTGTAATTATCTTTTCAAACTAAATATTAC　　

VCF

chr1	109817590	[known_CEBP_binding_increase]	G	T
chr10	23508363	[known_FOXA2_binding_decrease]	A	G
chr16	52599188	[known_FOXA1_binding_increase]	C	T
chr16	209709	[known_GATA1_binding_increase]	T	C　　

BED

chr1	109817090	109818090	.	0	*
chr10	23507863	23508863	.	0	*
chr16	52598688	52599688	.	0	*
chr16	209209	210209	.	0	*

 

網頁結果：

 

由於motif的binding是一個機率學的東西，因此理論上單個鹼基的改變會致使TF在該motif上binding的增長和減小，網頁裏面就是想預測出來某個特定的variant會致使哪些TF的binding出現顯著差別。

可是仍是不懂圖裏面的chromatin feature和log2FC是什麼意思！

如何理解這個Functional significance score，以及整個table？　　

 

關於模型的訓練：

數據

ENCODE and Roadmap Epigenomics data were used for labeling and the HG19 human genome was used for input sequences. The data is splitted to training, validation and test sets. The genomic regions are splitted to 200bp bins and labeled according to chromatin profiles. We kept the bins that have at least one TF binding event (note that TF binding event is measured by any overlap with a TF peak, not the >50% overlap criterion used for labeling).

DEPENDENCIES

1. Install CUDA driver. A high-end NVIDIA CUDA compatible graphics card with enough memory is required to train the model. I use Tesla K20m with 5Gb memory for training the model.

2. Installing torch and basic package dependencies following instructions from
http://torch.ch/docs/getting-started.html
You may need to install cmake first if you do not have it installed. It is highly recommended to link against OpenBLAS or other optimized BLAS library when building torch, since it makes a huge difference in performance while running on CPU.

3. Install torch packages required for training only: cutorch, cunn, mattorch. You may install through `luarock install [PACKAGE_NAME]` command. Note mattorch requires matlab. If you do not have matlab, you may try out https://github.com/soumith/matio-ffi.torch and change 1_data.lua to use matio instead (IMPORTANT: if you use matio, place remove the ":tr\
anspose(3,1)" and "transpose(2,1)" operation in 1_data.lua. The dimesions have been correctly handled by matio.).

這種必需要用顯卡了，N卡，用的CUDA平臺。二手的Tesla K20m with 5Gb memory大概兩千多。

Usage Example

th main.lua -save results

The output folder will be under ./results . The folder will inlcude the model file as well as log files for monitoring training progress.

You can specify various parameters for main.lua e.g. set learning rate by -LearningRate. Take a look at main.lua for the options.

Short explanation of the code: 1_data.lua reads the training and validation data; 2_model.lua specify the model; 3_loss.lua specify the loss function; 4_train.lua do the training.

能夠看到訓練model的任務是用lua語言來寫的。

 

關於Sequence Profiler功能

"in silico saturated mutagenesis" analysis for discovering informative sequence features within any sequence. 

預測所提供序列上每個位點發生各類突變時的log2FC的變化狀況。

 

 

參考：

深度學習在生物信息學上的應用

zhoujian最新的文章：Nature Genetics | 深度學習揭示大量自閉症相關非編碼區突變 

Whole-genome deep-learning analysis identifies contribution of noncoding mutations to autism risk

DeepSignal: 利用深度學習從 Nanopore 測序 reads 中檢測 DNA 甲基化狀態

DeepSignal: detecting DNA methylation state from Nanopore sequencing reads using deep-learning

第168期 | DEEPCPG：利用深度學習精確預測單細胞DNA甲基化狀態

DeepCpG: accurate prediction of single-cell DNA methylation states using deep learning