主要參考自：Hail | GWAS Tutorial
本筆記本旨在提供Hail功能的概述，重點(diǎn)是操作和查詢遺傳數(shù)據(jù)集的功能。我們進(jìn)行了全基因組SNP關(guān)聯(lián)測(cè)試，并證明了需要控制由群體分層引起的混雜。
1、安裝和環(huán)境準(zhǔn)備
#安裝 Jupyter
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh
source miniconda3/bin/activate
conda create -n python36 python=3.6
conda activate python36
pip install hail==0.2.11
pip install jupyterlab
sudo apt-get install -y openjdk-8-jre-headless g++ libopenblas-base liblapack3
#配置 Jupyter
touch venv/conf.py
#創(chuàng)建密碼
PASSWD=$(python -c 'from notebook.auth import passwd; print(passwd("jup"))')
echo "c.NotebookApp.password = u'${PASSWD}'"
vi conf.py
c.NotebookApp.open_browser = False
c.NotebookApp.port = 8818
c.NotebookApp.password = u'sha1:｛ sha1 值｝' # ${PASSWD} 替換為實(shí)際的 sha1 值
jupyter lab --config ./conf.py
# 新建一個(gè)python筆記本，初始化
import hail as hl
hl.init()
# ############
using hail jar at /home/ubuntu/miniconda3/envs/python36/lib/python3.6/site-packages/hail/hail-all-spark.jar
Running on Apache Spark version 2.2.3
SparkUI available at http://10.0.8.13:4040
Welcome to
__ __ <>__
/ /_/ /__ __/ /
/ __ / _ `/ / /
/_/ /_/\_,_/_/_/ version 0.2.11-cf54f08305d1
LOGGING: writing to /home/ubuntu/hail/hail-20220203-1927-0.2.11-cf54f08305d1.log
根據(jù)提示，看到的spark大數(shù)據(jù)的網(wǎng)頁(yè)管理界面
在使用 Hail 之前，我們導(dǎo)入了一些標(biāo)準(zhǔn)的 Python 庫(kù)，以便在整個(gè)筆記本中使用。
from hail.plot import show
from pprint import pprint
hl.plot.output_notebook()
這是一個(gè)類似ggplot這種包的python繪圖包。
# 下載1kg數(shù)據(jù)，好像還是谷歌的網(wǎng)址，竟然下載成功了
# 我們使用公共1000基因組數(shù)據(jù)集的一小部分，該數(shù)據(jù)集是通過(guò)將完整VCF中的基因分型SNP縮減到約20 MB采樣而創(chuàng)建的。我們還將集成來(lái)自單獨(dú)文本文件的示例和變體元數(shù)據(jù)。
# 這些文件由Hail團(tuán)隊(duì)托管在公共Google存儲(chǔ)桶中;以下命令下載該數(shù)據(jù)到本地。
hl.utils.get_1kg('data/')
# ####################
2022-02-03 19:59:22 Hail: INFO: downloading 1KG VCF ...
Source: https://storage.googleapis.com/hail-tutorial/1kg.vcf.bgz
2022-02-03 19:59:24 Hail: INFO: importing VCF and writing to matrix table...
2022-02-03 19:59:29 Hail: INFO: Coerced sorted dataset
2022-02-03 19:59:40 Hail: INFO: wrote matrix table with 10879 rows and 284 columns in 16 partitions to data/1kg.mt
2022-02-03 19:59:40 Hail: INFO: downloading 1KG annotations ...
Source: https://storage.googleapis.com/hail-tutorial/1kg_annotations.txt
2022-02-03 19:59:40 Hail: INFO: downloading Ensembl gene annotations ...
Source: https://storage.googleapis.com/hail-tutorial/ensembl_gene_annotations.txt
2022-02-03 19:59:42 Hail: INFO: Done!
# 查看文件情況
!ls -lh data
# total 18M
-rw-r--r-- 1 ubuntu ubuntu 100K Feb 3 19:59 1kg_annotations.txt
drwxrwxr-x 7 ubuntu ubuntu 4.0K Feb 3 19:59 1kg.mt
-rw-r--r-- 1 ubuntu ubuntu 16M Feb 3 19:59 1kg.vcf.bgz
-rw-r--r-- 1 ubuntu ubuntu 2.5M Feb 3 19:59 ensembl_gene_annotations.txt
2、先練練手，熟悉操作
環(huán)境配置好，可以愉快地學(xué)習(xí)啦！
從 VCF 導(dǎo)入數(shù)據(jù)
VCF 文件中的數(shù)據(jù)自然表示為 Hail MatrixTable。通過(guò)首先導(dǎo)入VCF文件，然后以Hail的文件格式寫(xiě)入生成的 MatrixTable，這樣對(duì)VCF數(shù)據(jù)的所有下游操作都將快得多。
hl.import_vcf('data/1kg.vcf.bgz').write('data/1kg.mt', overwrite=True)
# 2022-02-03 20:09:22 Hail: INFO: Coerced sorted dataset
2022-02-03 20:09:28 Hail: INFO: wrote matrix table with 10879 rows and 284 columns in 1 partition to data/1kg.mt
接下來(lái)，我們讀取寫(xiě)入的文件，分配變量mt（matrix table）。
mt = hl.read_matrix_table('data/1kg.mt')
了解我們的數(shù)據(jù)
重要的是要有簡(jiǎn)單的方法來(lái)切片、切塊、查詢和匯總數(shù)據(jù)集。下面演示了其中一些功能：
rows 方法可用于獲取包含 MatrixTable 中所有行字段的表。
我們可以與選擇一起使用以拉出5個(gè)變異。該方法采用引用表中字段名稱的字符串或 Hail Expression。在這里，我們將參數(shù)留空，以僅保留行鍵字段和 。rows select locus alleles
mt.rows().select().show(5)
# 或者 mt.row_key.show(5)
使用該方法顯示多屬性。show
mt.s.show(5)
要查看前幾個(gè)基因型調(diào)用，我們可以將entries與
select和take一起使用。take方法將前 n 行收集到一個(gè)列表中?；蛘撸覀兛梢允褂?code>show方法以表格格式將前n行打印到控制臺(tái)。
嘗試在下面的單元格中take更改為show。
mt.entry.take(5)
# ######
[Struct(GT=Call(alleles=[0, 0], phased=False), AD=[4, 0], DP=4, GQ=12, PL=[0, 12, 147]),
Struct(GT=Call(alleles=[0, 0], phased=False), AD=[8, 0], DP=8, GQ=24, PL=[0, 24, 315]),
Struct(GT=Call(alleles=[0, 0], phased=False), AD=[8, 0], DP=8, GQ=23, PL=[0, 23, 230]),
Struct(GT=Call(alleles=[0, 0], phased=False), AD=[7, 0], DP=7, GQ=21, PL=[0, 21, 270]),
Struct(GT=Call(alleles=[0, 0], phased=False), AD=[5, 0], DP=5, GQ=15, PL=[0, 15, 205])]
改為show看起來(lái)表格更好看：
添加列字段
Hail MatrixTable 可以具有任意數(shù)量的行字段和列字段，用于存儲(chǔ)與每行和每列關(guān)聯(lián)的數(shù)據(jù)。注釋通常是任何遺傳研究的關(guān)鍵部分。列字段用于存儲(chǔ)有關(guān)樣本表型、祖先、性別和協(xié)變量的信息。行字段可用于存儲(chǔ)成員基因和功能影響等信息，以便在QC或分析中使用。
在本教程中，我們將演示如何獲取文本文件并使用它來(lái)注釋 MatrixTable 中的列。
提供的文件包含樣本 ID、人口（國(guó)家）和"人口（地域）"名稱、樣本性別以及兩種模擬表型（二分類，或離散）。
此文件可以通過(guò)import_table導(dǎo)入到 Hail 中。此函數(shù)生成一個(gè) Table 對(duì)象。可以將其視為不受計(jì)算機(jī)上內(nèi)存限制的Pandas或R數(shù)據(jù)幀 - 在幕后，它用Spark。
# 導(dǎo)入
table = (hl.import_table('data/1kg_annotations.txt', impute=True)
.key_by('Sample'))
# #############
2022-02-07 15:09:18 Hail: INFO: Reading table to impute column types
2022-02-07 15:09:18 Hail: INFO: Finished type imputation
Loading column 'Sample' as type 'str' (imputed)
Loading column 'Population' as type 'str' (imputed)
Loading column 'SuperPopulation' as type 'str' (imputed)
Loading column 'isFemale' as type 'bool' (imputed)
Loading column 'PurpleHair' as type 'bool' (imputed)
Loading column 'CaffeineConsumption' as type 'int32' (imputed)
審視表格結(jié)構(gòu)的一個(gè)好方法是查看其schema。
table.describe()
# #############
----------------------------------------
Global fields:
None
----------------------------------------
Row fields:
'Sample': str
'Population': str
'SuperPopulation': str
'isFemale': bool
'PurpleHair': bool
'CaffeineConsumption': int32
----------------------------------------
Key: ['Sample']
----------------------------------------
若要查看前幾個(gè)值，請(qǐng)使用方法show:
table.show(width=100) # 100個(gè)字符？
# ########
現(xiàn)在，我們將使用此表將示例批注添加到數(shù)據(jù)集中，并將批注存儲(chǔ)在 MatrixTable 的列字段中。首先，我們將打印現(xiàn)有的列架構(gòu)（類似R語(yǔ)言class?）：
print(mt.col.dtype)
# struct{s: str}
我們使用annotate_cols方法將表與包含數(shù)據(jù)集的 MatrixTable 聯(lián)接在一起。
mt = mt.annotate_cols(pheno = table[mt.s])
mt.col.describe()
# #############
--------------------------------------------------------
Type:
struct {
s: str,
pheno: struct {
Population: str,
SuperPopulation: str,
isFemale: bool,
PurpleHair: bool,
CaffeineConsumption: int32
}
}
--------------------------------------------------------
Source:
<hail.matrixtable.MatrixTable object at 0x7fa7f1a38278>
Index:
['column']
--------------------------------------------------------
查詢函數(shù)和Hail表達(dá)式語(yǔ)言
Hail 有許多有用的查詢函數(shù)，可用于收集數(shù)據(jù)集上的統(tǒng)計(jì)信息。這些查詢函數(shù)將** Hail 表達(dá)式作為參數(shù)**。
我們將首先查看表中信息的一些統(tǒng)計(jì)信息。aggregate方法可用于聚合表中的行。
counter是一個(gè)聚合函數(shù)，用于計(jì)算每個(gè)唯一元素的出現(xiàn)次數(shù)。我們可以使用它來(lái)看人口的分布，方法是為我們要計(jì)數(shù)的字段傳遞Hail表達(dá)式。
pprint(table.aggregate(hl.agg.counter(table.SuperPopulation)))
# ##########
{'AFR': 1018, 'AMR': 535, 'EAS': 617, 'EUR': 669, 'SAS': 661}
2022-02-07 15:22:59 Hail: INFO: Coerced sorted dataset
# 看起來(lái)這個(gè)數(shù)據(jù)集清理不完全，好在只是警告，沒(méi)報(bào)錯(cuò)。
stats是一個(gè)聚合函數(shù)，用于生成有關(guān)數(shù)值集合的一些有用的統(tǒng)計(jì)信息。我們可以用它來(lái)觀察咖啡因消費(fèi)表型的分布。
pprint(table.aggregate(hl.agg.stats(table.CaffeineConsumption)))
# 幾個(gè)統(tǒng)計(jì)數(shù)字
{'max': 10.0,
'mean': 3.983714285714278,
'min': -1.0,
'n': 3500,
'stdev': 1.7021055628070707,
但是，這些指標(biāo)并不能完全代表我們數(shù)據(jù)集中的樣本。原因如下：
table.count()
# 3500
mt.count_cols()
# 284
mt.count_cols()
# 10879
由于數(shù)據(jù)集中的樣本少于完整的千人基因組計(jì)劃中的樣本，因此我們需要查看數(shù)據(jù)集上的注釋。我們可以使用aggregate_cols來(lái)僅獲取數(shù)據(jù)集中示例的指標(biāo)。
mt.aggregate_cols(hl.agg.counter(mt.pheno.SuperPopulation))
# #########少了很多樣本的
{'AFR': 76, 'EAS': 72, 'AMR': 34, 'SAS': 55, 'EUR': 47}
pprint(mt.aggregate_cols(hl.agg.stats(mt.pheno.CaffeineConsumption)))
{'max': 9.0,
'mean': 4.415492957746476,
'min': 0.0,
'n': 284,
'stdev': 1.5777634274659174,
'sum': 1253.9999999999993}
最后幾個(gè)命令中演示的功能并不是什么特別新的東西：使用Pandas或R數(shù)據(jù)幀，甚至是Unix工具（如awk）來(lái)解決這些問(wèn)題當(dāng)然不難。但是 Hail 可以使用相同的接口和查詢語(yǔ)言來(lái)分析要大得多的集，例如變異集。
在這里，我們計(jì)算 12 個(gè)可能的唯一 SNP 中每個(gè) SNP 的計(jì)數(shù)（參考?jí)A基有 4 個(gè) * 變異堿基 3 個(gè)=12個(gè)選項(xiàng)）。
為此，我們需要獲取每個(gè)變體的變異等位基因，然后計(jì)算每個(gè)唯一ref / alt對(duì)的出現(xiàn)次數(shù)。這可以通過(guò)Hail的counter功能來(lái)完成。
snp_counts = mt.aggregate_rows(hl.agg.counter(hl.Struct(ref=mt.alleles[0], alt=mt.alleles[1])))
pprint(snp_counts)
# ##########
{Struct(ref='A', alt='C'): 451,
Struct(ref='C', alt='A'): 494,
Struct(ref='C', alt='G'): 150,
Struct(ref='C', alt='T'): 2418,
Struct(ref='G', alt='A'): 2367,
Struct(ref='A', alt='T'): 75,
Struct(ref='T', alt='G'): 466,
Struct(ref='G', alt='C'): 111,
Struct(ref='G', alt='T'): 477,
Struct(ref='A', alt='G'): 1929,
Struct(ref='T', alt='A'): 77,
Struct(ref='T', alt='C'): 1864}
可以使用Python的Counter類按降序列出計(jì)數(shù)。
from collections import Counter
counts = Counter(snp_counts)
counts.most_common()
# ############
[(Struct(ref='C', alt='T'), 2418),
(Struct(ref='G', alt='A'), 2367),
(Struct(ref='A', alt='G'), 1929),
(Struct(ref='T', alt='C'), 1864),
(Struct(ref='C', alt='A'), 494),
(Struct(ref='G', alt='T'), 477),
(Struct(ref='T', alt='G'), 466),
(Struct(ref='A', alt='C'), 451),
(Struct(ref='C', alt='G'), 150),
(Struct(ref='G', alt='C'), 111),
(Struct(ref='T', alt='A'), 77),
(Struct(ref='A', alt='T'), 75)]
很高興看到我們實(shí)際上可以從這個(gè)小數(shù)據(jù)集中發(fā)現(xiàn)一些生物學(xué)的東西：我們看到這些頻率成對(duì)出現(xiàn)。C / T和G / A實(shí)際上是相同的突變，只是從相反的鏈中觀察。同樣，T/A 和 A/T 在相反的鏈上是相同的突變。C/T 和 A/T SNP 的頻率之間有 30 倍的差異。為什么？
相同的Python，R和Unix工具也可以完成這項(xiàng)工作，但我們開(kāi)始碰壁 - 最新的gnomaD版本發(fā)布了大約2.5億個(gè)變體，并且無(wú)法在一臺(tái)計(jì)算機(jī)上內(nèi)存中。
基因型呢？Hail可以查詢數(shù)據(jù)集中所有基因型的集合，即使對(duì)于我們這個(gè)小數(shù)據(jù)集來(lái)說(shuō)，這也越來(lái)越大。我們的 284 個(gè)樣本和 10，000 個(gè)變異可產(chǎn)生 1000 萬(wàn)個(gè)獨(dú)特的基因型。gnomAD數(shù)據(jù)集有大約5萬(wàn)億個(gè)獨(dú)特的基因型。
Hail繪制函數(shù)允許Hail字段作為參數(shù)，因此我們可以直接在此處傳遞 DP 字段。如果未設(shè)置范圍和條柱參數(shù)，則此函數(shù)將根據(jù)字段的最小值和最大值計(jì)算范圍，并使用默認(rèn)的 50 個(gè)柱子。
p = hl.plot.histogram(mt.DP, range=(0,30), bins=30, title='DP Histogram', legend='DP')
show(p)
3、質(zhì)控
分析師將大部分時(shí)間花在測(cè)序數(shù)據(jù)集的QC上。QC是一個(gè)迭代過(guò)程，每個(gè)項(xiàng)目都是不同的：QC沒(méi)有"按鈕"解決方案。每次Broad收集一組新的樣本時(shí)，它都會(huì)發(fā)現(xiàn)新的批次效應(yīng)。但是，通過(guò)實(shí)踐開(kāi)放科學(xué)并與他人討論QC流程和決策，我們可以作為一個(gè)社區(qū)建立一套最佳實(shí)踐。
QC完全基于理解數(shù)據(jù)集屬性的能力。Hail 試圖通過(guò)提供sample_qc函數(shù)來(lái)簡(jiǎn)化此操作，該函數(shù)生成一組有用的指標(biāo)并將其存儲(chǔ)在列字段中。
mt.col.describe()
# ##########
--------------------------------------------------------
Type:
struct {
s: str,
pheno: struct {
Population: str,
SuperPopulation: str,
isFemale: bool,
PurpleHair: bool,
CaffeineConsumption: int32
}
}
--------------------------------------------------------
Source:
<hail.matrixtable.MatrixTable object at 0x7fa7f1a38278>
Index:
['column']
--------------------------------------------------------
mt = hl.sample_qc(mt)
--------------------------------------------------------
Type:
struct {
s: str,
pheno: struct {
Population: str,
SuperPopulation: str,
isFemale: bool,
PurpleHair: bool,
CaffeineConsumption: int32
},
sample_qc: struct {
dp_stats: struct {
mean: float64,
stdev: float64,
min: float64,
max: float64
},
gq_stats: struct {
mean: float64,
stdev: float64,
min: float64,
max: float64
},
call_rate: float64,
n_called: int64,
n_not_called: int64,
n_hom_ref: int64,
n_het: int64,
n_hom_var: int64,
n_non_ref: int64,
n_singleton: int64,
n_snp: int64,
n_insertion: int64,
n_deletion: int64,
n_transition: int64,
n_transversion: int64,
n_star: int64,
r_ti_tv: float64,
r_het_hom_var: float64,
r_insertion_deletion: float64
}
}
mt.col.describe()
--------------------------------------------------------
Source:
<hail.matrixtable.MatrixTable object at 0x7fa7f09bd240>
Index:
['column']
--------------------------------------------------------
繪制 QC 指標(biāo)是一個(gè)很好的起點(diǎn)。
p = hl.plot.histogram(mt.sample_qc.call_rate, range=(.88,1), legend='Call Rate')
show(p)
p = hl.plot.histogram(mt.sample_qc.gq_stats.mean, range=(10,70), legend='Mean Sample GQ')
show(p)
通常，這些指標(biāo)是相關(guān)的。
p = hl.plot.scatter(mt.sample_qc.dp_stats.mean, mt.sample_qc.call_rate, xlabel='Mean DP', ylabel='Call Rate')
show(p)
從數(shù)據(jù)集中刪除異常值通常會(huì)改善關(guān)聯(lián)結(jié)果暫時(shí)沒(méi)搞清楚DP，AD是啥意思，繼續(xù)。我們可以進(jìn)行"任意(這里應(yīng)該是嘗試多次的意思，歪國(guó)人喜歡搞笑)"截止，并使用它們來(lái)過(guò)濾：
mt = mt.filter_cols((mt.sample_qc.dp_stats.mean >= 4) & (mt.sample_qc.call_rate >= 0.97))
print('After filter, %d/284 samples remain.' % mt.count_cols())
# After filter, 250/284 samples remain.
接下來(lái)是基因型QC。在讀數(shù)不在應(yīng)該的地方過(guò)濾掉基因型是個(gè)好主意：如果我們找到一個(gè)稱為純合子參考的基因型>10%的alt reads，或者稱為雜合子的基因型，沒(méi)有接近1：1的ref / alt平衡，則很可能是一個(gè)錯(cuò)誤。
在像1KG這樣的低深度數(shù)據(jù)集中(PS..有點(diǎn)改變認(rèn)知，1kG不是高深度度測(cè)序的，可是全人類的研究標(biāo)準(zhǔn)呀)，很難使用此指標(biāo)檢測(cè)不良基因型，因?yàn)?與10的reads比率可以很容易地用二項(xiàng)采樣來(lái)解釋。但是，在深度數(shù)據(jù)集中，10：100 的讀取比率肯定會(huì)引起關(guān)注！
ab = mt.AD[1] / hl.sum(mt.AD)
filter_condition_ab = ((mt.GT.is_hom_ref() & (ab <= 0.1)) |
(mt.GT.is_het() & (ab >= 0.25) & (ab <= 0.75)) |
(mt.GT.is_hom_var() & (ab >= 0.9)))
fraction_filtered = mt.aggregate_entries(hl.agg.fraction(~filter_condition_ab))
print(f'Filtering {fraction_filtered * 100:.2f}% entries out of downstream analysis.')
mt = mt.filter_entries(filter_condition_ab)
# Filtering 3.63% entries out of downstream analysis.
變異QC 更是一樣的：我們可以使用 variant_qc 函數(shù)來(lái)生成各種有用的統(tǒng)計(jì)數(shù)據(jù)，繪制它們并進(jìn)行篩選。
mt = hl.variant_qc(mt)
mt.row.describe()
--------------------------------------------------------
Type:
struct {
locus: locus<GRCh37>,
alleles: array<str>,
rsid: str,
qual: float64,
filters: set<str>,
info: struct {
AC: array<int32>,
AF: array<float64>,
AN: int32,
BaseQRankSum: float64,
ClippingRankSum: float64,
DP: int32,
DS: bool,
FS: float64,
HaplotypeScore: float64,
InbreedingCoeff: float64,
MLEAC: array<int32>,
MLEAF: array<float64>,
MQ: float64,
MQ0: int32,
MQRankSum: float64,
QD: float64,
ReadPosRankSum: float64,
set: str
},
variant_qc: struct {
AC: array<int32>,
AF: array<float64>,
AN: int32,
homozygote_count: array<int32>,
dp_stats: struct {
mean: float64,
stdev: float64,
min: float64,
max: float64
},
gq_stats: struct {
mean: float64,
stdev: float64,
min: float64,
max: float64
},
n_called: int64,
n_not_called: int64,
call_rate: float32,
n_het: int64,
n_non_ref: int64,
het_freq_hwe: float64,
p_value_hwe: float64
}
}
--------------------------------------------------------
Source:
<hail.matrixtable.MatrixTable object at 0x7fa7f09dd6a0>
Index:
['row']
--------------------------------------------------------
這些統(tǒng)計(jì)數(shù)據(jù)實(shí)際上看起來(lái)相當(dāng)不錯(cuò)：我們不需要過(guò)濾這個(gè)數(shù)據(jù)集。但是，大多數(shù)數(shù)據(jù)集都需要經(jīng)過(guò)深思熟慮的質(zhì)量控制。filter_rows方法可以提供幫助！
4、開(kāi)始GWAS！
首先，我們需要限制為以下變體：
常見(jiàn)（我們將使用 **1% **的截止值）
距離哈代-溫伯格均衡不遠(yuǎn)，來(lái)發(fā)現(xiàn)測(cè)序錯(cuò)誤
mt = mt.filter_rows(mt.variant_qc.AF[1] > 0.01)
mt = mt.filter_rows(mt.variant_qc.p_value_hwe > 1e-6)
print('Samples: %d Variants: %d' % (mt.count_cols(), mt.count_rows()))
# Samples: 250 Variants: 7779
這些過(guò)濾器刪除了大約15%的位點(diǎn)（我們從10，000多個(gè)位點(diǎn)開(kāi)始）。這并不代表大多數(shù)測(cè)序數(shù)據(jù)集！我們已經(jīng)對(duì)整整一千個(gè)基因組數(shù)據(jù)集進(jìn)行了縮減采樣，以包括比我們偶然預(yù)期的更常見(jiàn)的變體。
在 Hail 中，關(guān)聯(lián)檢驗(yàn)接受樣本表型和協(xié)變量的列字段。我們已經(jīng)在數(shù)據(jù)集中獲得了感興趣的表型（咖啡因消耗）：
gwas = hl.linear_regression_rows(y=mt.pheno.CaffeineConsumption,
x=mt.GT.n_alt_alleles(),
covariates=[1.0])
gwas.row.describe()
--------------------------------------------------------
Type:
struct {
locus: locus<GRCh37>,
alleles: array<str>,
n: int32,
sum_x: float64,
y_transpose_x: float64,
beta: float64,
standard_error: float64,
t_stat: float64,
p_value: float64
}
--------------------------------------------------------
Source:
<hail.table.Table object at 0x7fa7f0945828>
Index:
['row']
--------------------------------------------------------
查看上述打印輸出的底部，您可以看到線性回歸為 beta、標(biāo)準(zhǔn)誤差、t 統(tǒng)計(jì)量和 p 值添加了新的行字段。
Hail使結(jié)果可視化變得容易！讓我們做一個(gè)曼哈頓圖：
這看起來(lái)不像一個(gè)天際線。讓我們檢查一下我們的GWAS是否使用Q-Q（分位數(shù)-分位數(shù)）圖得到了很好的控制。
p = hl.plot.qq(gwas.p_value)
show(p)
5、混淆！
觀察到的 p 值會(huì)全部偏離預(yù)期。要么我們數(shù)據(jù)集中的每個(gè)SNP都與咖啡因的攝入有因果關(guān)系（不太可能），要么有一個(gè)混雜因素。
我們沒(méi)有告訴你，但樣本祖先實(shí)際上被用來(lái)模擬這種表型。這導(dǎo)致表型的分層分布。解決方案是在我們的回歸中將祖先作為協(xié)變量包括在內(nèi)。
linear_regression_rows函數(shù)還可以采用列字段作為協(xié)變量使用。我們已經(jīng)用報(bào)告的祖先注釋了我們的樣本，但由于人為錯(cuò)誤，對(duì)這些標(biāo)簽持懷疑態(tài)度是件好事。基因組沒(méi)有這個(gè)問(wèn)題！我們將通過(guò)使用報(bào)告的祖先，而是通過(guò)在我們的模型中包含計(jì)算的主成分來(lái)作為遺傳祖先。
pca 函數(shù)將特征值生成為列表，將示例 PC 生成為表，并且還可以在詢問(wèn)時(shí)生成變體加載。hwe_normalized_pca函數(shù)也這樣做，使用HWE歸一化基因型進(jìn)行PCA。
eigenvalues, pcs, _ = hl.hwe_normalized_pca(mt.GT)
# ##################
2022-02-08 09:34:10 Hail: INFO: hwe_normalized_pca: running PCA using 7771 variants.
2022-02-08 09:34:15 Hail: INFO: pca: running PCA with 10 components...
pprint(eigenvalues)
# ##########
[18.078002475659527,
9.978057054208994,
3.536048746504391,
2.6583611350017096,
1.5966283377451316,
1.5408505989074688,
1.5064316096906518,
1.4747771722760974,
1.4667919317402691,
1.4458363475579015]
pcs.show(5, width=100)
現(xiàn)在我們已經(jīng)為每個(gè)樣本提供了主成分，我們不妨繪制它們！人類歷史在遺傳數(shù)據(jù)集中產(chǎn)生了強(qiáng)烈的影響。即使使用50MB的測(cè)序數(shù)據(jù)集，我們也可以恢復(fù)主要人群。
mt = mt.annotate_cols(scores = pcs[mt.s].scores)
p = hl.plot.scatter(mt.scores[0],
mt.scores[1],
label=mt.pheno.SuperPopulation,
title='PCA', xlabel='PC1', ylabel='PC2')
show(p)
現(xiàn)在，我們可以重新運(yùn)行線性回歸，控制樣本性別和前幾個(gè)主成分。我們將像以前一樣使用輸入變量替代等位基因的數(shù)量來(lái)執(zhí)行此操作，并再次使用輸入變量從PL字段導(dǎo)出的基因型劑量。
gwas = hl.linear_regression_rows(
y=mt.pheno.CaffeineConsumption,
x=mt.GT.n_alt_alleles(),
covariates=[1.0, mt.pheno.isFemale, mt.scores[0], mt.scores[1], mt.scores[2]])
# 2022-02-08 09:49:57 Hail: INFO: linear_regression_rows: running on 250 samples for 1 response variable y,
with input variable x, and 5 additional covariates...
我們將首先制作一個(gè)Q-Q圖來(lái)評(píng)估膨脹...
這種形狀表明一項(xiàng)控制良好（但不是特別有效）的研究?，F(xiàn)在是曼哈頓的情節(jié)：
p = hl.plot.manhattan(gwas.p_value)
show(p)
我們發(fā)現(xiàn)了一個(gè)咖啡因的消費(fèi)點(diǎn)！現(xiàn)在只需應(yīng)用Hail的Nature論文功能即可發(fā)布結(jié)果。
開(kāi)個(gè)玩笑，這個(gè)功能要到Hail 1.0才會(huì)登陸！
罕見(jiàn)變異分析
在這里，我們將演示如何使用表達(dá)式語(yǔ)言按行和列字段中的任何任意屬性進(jìn)行分組和計(jì)數(shù)。Hail 還實(shí)現(xiàn)了序列核心關(guān)聯(lián)測(cè)檢驗(yàn)（SKAT）。
entries = mt.entries()
results = (entries.group_by(pop = entries.pheno.SuperPopulation, chromosome = entries.locus.contig)
.aggregate(n_het = hl.agg.count_where(entries.GT.is_het())))
results.show()
我們使用 MatrixTable.entries 方法將矩陣表轉(zhuǎn)換為表（每個(gè)變體的每個(gè)樣本都有一行）。在這種表示中，很容易對(duì)我們喜歡的任何字段進(jìn)行聚合，這通常是罕見(jiàn)變體分析的第一步。
如果我們想按次要等位基因頻率箱和頭發(fā)顏色分組，并計(jì)算平均GQ，該怎么辦？
我們已經(jīng)證明，通過(guò)幾個(gè)任意統(tǒng)計(jì)數(shù)據(jù)很容易聚合。此特定示例可能無(wú)法提供特別有用的信息，但相同的模式可用于檢測(cè)罕見(jiàn)變異的影響：
entries = entries.annotate(maf_bin = hl.if_else(entries.info.AF[0]<0.01, "< 1%",
hl.if_else(entries.info.AF[0]<0.05, "1%-5%", ">5%")))
results2 = (entries.group_by(af_bin = entries.maf_bin, purple_hair = entries.pheno.PurpleHair)
.aggregate(mean_gq = hl.agg.stats(entries.GQ).mean,
mean_dp = hl.agg.stats(entries.DP).mean))
results2.show()
# module 'hail' has no attribute 'if_else'
# 有點(diǎn)奇怪的報(bào)錯(cuò)，暫未找到原因
按功能類別（同義詞、錯(cuò)義或功能喪失）計(jì)算每個(gè)基因的雜合子基因型數(shù)量，以估計(jì)每個(gè)基因的功能約束
計(jì)算病例中每個(gè)基因的單例功能喪失突變數(shù)，并進(jìn)行對(duì)照，以檢測(cè)與疾病有關(guān)的基因。
結(jié)語(yǔ)
恭喜！您已經(jīng)到了第一個(gè)教程的末尾。要了解有關(guān) Hail 的 API 和功能的更多信息，請(qǐng)查看其他教程。您可以查看 Python API 以獲取有關(guān)其他 Hail 函數(shù)的文檔。如果您將Hail用于自己的科學(xué)，我們很樂(lè)意在Zulip聊天或論壇上收到您的來(lái)信。
作為參考，下面是合并到一個(gè)腳本中的所有教程終結(jié)點(diǎn)的完整流程。
table = hl.import_table('data/1kg_annotations.txt', impute=True).key_by('Sample')
mt = hl.read_matrix_table('data/1kg.mt')
mt = mt.annotate_cols(pheno = table[mt.s])
mt = hl.sample_qc(mt)
mt = mt.filter_cols((mt.sample_qc.dp_stats.mean >= 4) & (mt.sample_qc.call_rate >= 0.97))
ab = mt.AD[1] / hl.sum(mt.AD)
filter_condition_ab = ((mt.GT.is_hom_ref() & (ab <= 0.1)) |
(mt.GT.is_het() & (ab >= 0.25) & (ab <= 0.75)) |
(mt.GT.is_hom_var() & (ab >= 0.9)))
mt = mt.filter_entries(filter_condition_ab)
mt = hl.variant_qc(mt)
mt = mt.filter_rows(mt.variant_qc.AF[1] > 0.01)
eigenvalues, pcs, _ = hl.hwe_normalized_pca(mt.GT)
mt = mt.annotate_cols(scores = pcs[mt.s].scores)
gwas = hl.linear_regression_rows(
y=mt.pheno.CaffeineConsumption,
x=mt.GT.n_alt_alleles(),
covariates=[1.0, mt.pheno.isFemale, mt.scores[0], mt.scores[1], mt.scores[2]])
