Spark 高级分析:第十章第5节 示例：从1000个基因组项目中查询基因型

问题导读

1.如何获取示例数据？
2.要将示例中的数据做哪些操作？如何操作？
3.基因组学中的许多计算是否很好地融入了Spark计算范例？




Spark 高级分析:第十章第4节 示例：从编码数据预测转录因子结合位点
http://www.aboutyun.com/forum.php?mod=viewthread&tid=26573



在这个例子中，我们将摄取完整的1000个基因组基因型数据集。首先，我们直接将原始数据下载到HDFS中，解压缩，然后运行ADAM作业将数据转换为Parquet格式。下面是一个示例命令；应该对所有染色体执行这个命令，并且可以在集群中并行化。
[mw_shl_code=bash,true]curl -s -L ftp://.../1000genomes/.../chr1.vcf.gz \
| gunzip \
| hadoop fs -put - /user/ds/genomics/1kg/vcf/chr1.vcf
export SPARK_JAR_PATH=hdfs:///path/to/spark.jar
adam/bin/adam-submit --conf spark.yarn.jar=$SPARK_JAR_PATH \
vcf2adam \
-coalesce 5 \
/user/ds/genomics/1kg/vcf/chr1.vcf \
/user/ds/genomics/1kg/parquet/chr1[/mw_shl_code]
请注意我们如何指定-coalesce 5；这将确保映射任务将数据压缩为少量的大Parquet文件。然后，从ADAM shell中加载并检查一个对象，例如：
[mw_shl_code=scala,true]import org.bdgenomics.adam.rdd.ADAMContext._
import org.bdgenomics.formats.avro.Genotype
val genotypesRDD = sc.adamLoad[Genotype, Nothing](
"/user/ds/genomics/1kg/parquet")
val gt = genotypesRDD.first()
...[/mw_shl_code]
假设我们想计算所有样本中每个重叠CTCF结合位点的变体基因组范围的次等位基因频率。我们必须将来自最后一部分的CTCF数据和来自1000基因组项目的基因型数据结合起来：
[mw_shl_code=scala,true]val ctcfRDD = sc.adamNarrowPeakFeatureLoad(
"/user/ds/genomics/chip-seq/GM12878.ChIP-seq.CTCF.narrowPeak")
val filtered = (BroadcastRegionJoin.partitionAndJoin(
sc, ctcfRDD, genotypesRDD)
.map(_._2))[/mw_shl_code]
我们还需要一个功能，将获取Genotype类型参数并计算计数的参考/替代等位基因：
[mw_shl_code=scala,true]def genotypeToAlleleCounts(gt: Genotype): (Variant, (Int, Int)) = {
val counts = gt.getAlleles.map(allele match {
case GenotypeAllele.Ref => (1, 0)
case GenotypeAllele.Alt => (0, 1)
case _ => (0, 0)
}).reduce((x, y) => (x._1 + y._1, x._2 + y._2))
(gt.getVariant, (counts._1, counts._2))
}[/mw_shl_code]
最后，我们生成RDD[(Variant, (Int, Int))]并执行聚合：
[mw_shl_code=scala,true]val counts = filtered.map(genotypeToAlleleCounts)
val countsByVariant = counts.reduceByKey(
(x, y) => (x._1 + y._1, x._2 + y._2))
val mafByVariant = countsByVariant.map(tup => {
val (v, (r, a)) = tup
val n = r + a
(v, math.min(r, a).toDouble / n)
})[/mw_shl_code]
遍历整个数据集是一个相当大的操作。因为我们只从基因型数据中访问少数字段，所以它肯定会受益于谓词下推和投影，我们将其作为练习留给读者。


下一步
基因组学中的许多计算很好地融入了Spark计算范例。在执行临时分析时，像ADAM这样的项目提供的最有价值的贡献是表示底层分析对象（以及转换工具）的Avro模式集。我们看到，一旦数据被转换成相应的Avro模式，许多大规模计算就变得相对容易表达和分布。

虽然可能仍然缺乏对Hadoop/Spark进行科学研究的工具，但是确实有一些项目可以帮助避免重新发明轮子。我们探讨了在ADAM中实现的核心功能，但是该项目已经为整个GATK最佳实践流水线提供了实现，包括BQSR、indel重新对齐和去重复。除了ADAM，许多机构已经与全球基因组学与健康联盟签约，该联盟已经开始产生他们自己的基因组分析模式。西奈山医学院的Hammerbacher实验室还开发了瓜卡莫尔（Guacamole），这是一套主要针对体细胞变异体的工具，呼吁癌症基因组学。所有这些工具都是具有自由Apache v2许可证的开放源码，因此如果你在自己的工作中开始使用它们，请考虑做出改进！

最新经典文章，欢迎关注公众号