about云开发

 找回密码
 立即注册

QQ登录

只需一步,快速开始

扫一扫,访问微社区

查看: 689|回复: 0

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

[复制链接]

101

主题

21

听众

3

收听

版主

Rank: 7Rank: 7Rank: 7

积分
1734
发表于 2019-1-11 08:54:39 | 显示全部楼层 |阅读模式
本帖最后由 feilong 于 2019-1-11 08:57 编辑

问题导读

1.
如何获取示例数据?
2.要将示例中的数据做哪些操作?如何操作?

3.
基因组学中的许多计算是否很好地融入了Spark计算范例





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



在这个例子中,我们将摄取完整的1000个基因组基因型数据集。首先,我们直接将原始数据下载到HDFS中,解压缩,然后运行ADAM作业将数据转换为Parquet格式。下面是一个示例命令;应该对所有染色体执行这个命令,并且可以在集群中并行化。
[Bash shell] 纯文本查看 复制代码
curl -s -L [url=ftp://.../1000genomes/.../chr1.vcf.gz]ftp://.../1000genomes/.../chr1.vcf.gz[/url] \
| 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

请注意我们如何指定-coalesce 5;这将确保映射任务将数据压缩为少量的大Parquet文件。然后,从ADAM shell中加载并检查一个对象,例如:
[Scala] 纯文本查看 复制代码
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()
...

假设我们想计算所有样本中每个重叠CTCF结合位点的变体基因组范围的次等位基因频率。我们必须将来自最后一部分的CTCF数据和来自1000基因组项目的基因型数据结合起来:
[Scala] 纯文本查看 复制代码
val ctcfRDD = sc.adamNarrowPeakFeatureLoad(
"/user/ds/genomics/chip-seq/GM12878.ChIP-seq.CTCF.narrowPeak")
val filtered = (BroadcastRegionJoin.partitionAndJoin(
sc, ctcfRDD, genotypesRDD)
.map(_._2))

我们还需要一个功能,将获取Genotype类型参数并计算计数的参考/替代等位基因:
[Scala] 纯文本查看 复制代码
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))
}

最后,我们生成RDD[(Variant, (Int, Int))]并执行聚合:
[Scala] 纯文本查看 复制代码
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)
})

遍历整个数据集是一个相当大的操作。因为我们只从基因型数据中访问少数字段,所以它肯定会受益于谓词下推和投影,我们将其作为练习留给读者。


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

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




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

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

关闭

推荐上一条 /3 下一条

QQ|小黑屋|about云开发-学问论坛|社区 ( 京ICP备12023829号

GMT+8, 2019-3-24 19:13 , Processed in 0.341858 second(s), 29 queries , Gzip On.

Powered by Discuz! X3.2 Licensed

快速回复 返回顶部 返回列表