A genomics processing engine and specialized file format built using Apache Avro, Apache Spark and Parquet. Apache 2 licensed.
Context: rust-bio/rust-htslib#40 ... I would like to have those as part of the output parquet file. Currently I'm trying to do that outside of ADAM as an extra step: https://github.com/brainstorm/htsget-indexer ... but I was wondering if this is a feature/option you had in mind at some point?
According to spec, those are calculated with some bit shifting... not sure if htsjdk has that in place, but it could/should come from there in the case of java, I guess.
htslib-rs has
.tell/.seek methods to retrieve those
Yeah, essentially the voffsets are to be paired with genomics coordinates, see here for a more straightforward and low level explanation: https://delagoya.github.io/rusty-bio/binary-file-parsing-in-bioinformatics
And htsjdk indeed seems to support those: https://github.com/samtools/htsjdk/blob/28dde96231dc9975a90ac2ca180c6a0991f2edcb/src/test/java/htsjdk/samtools/BAMSBIIndexerTest.java#L64
@heuermh Would it be acceptable to modify the
bdg.avdl to include a couple of fields such as coffset/uoffset? I would like to avoid chucking those into the additional attributes: https://github.com/bigdatagenomics/bdg-formats/blob/master/src/main/resources/avro/bdg.avdl#L446
I guess I'm on the fence, since those fields are specific to a flat file implementation that frankly we'd like to see go away. If an argument can be made that we need to keep those around in order to properly round-trip, and ADAM reads and writes those as a course of "normal business", then I might be convinced
Perhaps you might discuss the use case more? If alignment records are in Parquet on disk, and you need quick access to region based queries, we've persued those use cases with
IntervalRDD in Mango, and by writing Parquet partitioned by region, i.e. AlignmentDataset.saveAsPartitionedParquet(partitionSize: 100000)
Thanks Michael for the pointer and explanation. I do agree and thought about the inconvenience of putting those low level details (although just as an optional field)... the main use case would be https://samtools.github.io/hts-specs/htsget.html, where a particular genomic coordinate should return a JSON with matching reads and accompanying range of bytes from the original file.
I can see how Mango encompasses a lot of the fast indexing/interval tree ideas that htsget needs on the backend, so that's a big win/realization I did not have before... I'll explore a bit more Mango and see if there's a chance to expose those intervals as a htsget-like REST API.
With regard to the problem that China cannot be installed, it is suggested that Failed to execute goal on project adam-core-spark2_2.11: Could not resolve dependencies for project org.bdgenomics.adam:adam-core-spark2_2.11:jar:0.30.0-SNAPSHOT: The following artifacts could not be resolved: org.apache.avro:avro-mapred:jar:hadoop2:1.8.2, hopes to have a download source from China. If I give me the file, I can provide a compiled download source for adamundefineds mvn in China. For use by adam personnel in China to promote the spread of adam in China
New releases:
https://github.com/bigdatagenomics/adam version 0.30.0
https://github.com/bigdatagenomics/cannoli version 0.8.0
And related:
https://github.com/bigdatagenomics/bdg-formats version 0.15.0
https://github.com/bigdatagenomics/convert version 0.9.0
Any1 have an example of using adam to bin reads and run haplotypecaller over the bins? Wanting to run some avocado benchmarks against haplotypecaller, but i'd rather not use their spark implementation. Looking to just bin reads and run the single threaded version on each of the bins
val alignments = sc.loadAlignments("sample.bam")
alignments.saveAsPartitionedParquet("sample.partitioned.adam")
val bin = sc.loadPartitionedParquetAlignments("sample.partitioned.adam", regions = Seq(...))
Hmm I think that almost works. I think I would need to figure out how to "slop" those partitions in a way similar to this (https://bedtools.readthedocs.io/en/latest/content/tools/slop.html)
I sometimes get this error when running bwa and i'm not sure why:
Caused by: java.util.NoSuchElementException: key not found: sample
at scala.collection.MapLike$class.default(MapLike.scala:228)
at scala.collection.AbstractMap.default(Map.scala:59)
at scala.collection.MapLike$class.apply(MapLike.scala:141)
at scala.collection.AbstractMap.apply(Map.scala:59)
at org.bdgenomics.adam.models.ReadGroupDictionary.apply(ReadGroupDictionary.scala:133)
at org.bdgenomics.adam.converters.AlignmentConverter$$anonfun$convert$1$$anonfun$apply$1.apply(AlignmentConverter.scala:259)
at org.bdgenomics.adam.converters.AlignmentConverter$$anonfun$convert$1$$anonfun$apply$1.apply(AlignmentConverter.scala:257)It comes from this. Mid alignment the readgroupmap doesn't contain the sample name assigned.
Just a local test.
val readFile1 = "/home/james/genomic_data/reads/test_1.fq"
val readFile2 = "/home/james/genomic_data/reads/test_2.fq"
val ref = "/home/james/genomic_data/reference/human_g1k_v37.fasta"
val args = new BwaArgs()
args.indexPath = ref
args.sample = "sample"
val reads = sc.loadPairedFastqAsFragments(readFile1, readFile2)
val alignments = reads.alignWithBwa(args)
alignments.saveAsSam("/home/james/genomic_data/reads/test_1.bam", asSingleFile = true)
The cannoli bwa command line class does two extra steps after piping out to
https://github.com/bigdatagenomics/cannoli/blob/master/cli/src/main/scala/org/bdgenomics/cannoli/cli/Bwa.scala#L102
bwa, creating read groups and referenceshttps://github.com/bigdatagenomics/cannoli/blob/master/cli/src/main/scala/org/bdgenomics/cannoli/cli/Bwa.scala#L102
Interesting error
Exception in thread "main" java.lang.IllegalArgumentException: Fastq 1 (s3a://bgi-concordance-studies/run190704NovaSeq_R1_A501.fastq.gz.corrected.fastq.err_barcode_removed.fixed.fastq.gz) has 1023782058 reads, fastq 2 (s3a://bgi-concordance-studies/run190704NovaSeq_R2_A501.fastq.gz.corrected.fastq.err_barcode_removed.fixed.fastq.gz) has 1023782057 reads
at org.bdgenomics.adam.rdd.ADAMContext$$anonfun$loadPairedFastq$1.apply(ADAMContext.scala:2536)
at org.bdgenomics.adam.rdd.ADAMContext$$anonfun$loadPairedFastq$1.apply(ADAMContext.scala:2513)
at scala.Option.fold(Option.scala:158)
at org.apache.spark.rdd.Timer.time(Timer.scala:48)
at org.bdgenomics.adam.rdd.ADAMContext.loadPairedFastq(ADAMContext.scala:2513)
at org.bdgenomics.adam.rdd.ADAMContext$$anonfun$loadPairedFastqAsFragments$1.apply(ADAMContext.scala:3059)
at org.bdgenomics.adam.rdd.ADAMContext$$anonfun$loadPairedFastqAsFragments$1.apply(ADAMContext.scala:3059)
at scala.Option.fold(Option.scala:158)
at org.apache.spark.rdd.Timer.time(Timer.scala:48)
at org.bdgenomics.adam.rdd.ADAMContext.loadPairedFastqAsFragments(ADAMContext.scala:3057)
at org.myome.bigbioinf.Main$.main(Main.scala:39)
at org.myome.bigbioinf.Main.main(Main.scala)
at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
at java.lang.reflect.Method.invoke(Method.java:498)
at org.apache.spark.deploy.JavaMainApplication.start(SparkApplication.scala:52)
at org.apache.spark.deploy.SparkSubmit.org$apache$spark$deploy$SparkSubmit$$runMain(SparkSubmit.scala:845)
at org.apache.spark.deploy.SparkSubmit.doRunMain$1(SparkSubmit.scala:161)
at org.apache.spark.deploy.SparkSubmit.submit(SparkSubmit.scala:184)
at org.apache.spark.deploy.SparkSubmit.doSubmit(SparkSubmit.scala:86)
at org.apache.spark.deploy.SparkSubmit$$anon$2.doSubmit(SparkSubmit.scala:920)
at org.apache.spark.deploy.SparkSubmit$.main(SparkSubmit.scala:929)
at org.apache.spark.deploy.SparkSubmit.main(SparkSubmit.scala)The true number of reads is 1023782052
Somehow this is conditional is failing even though i have confirmed that the number of reads are the same
not sure if the bar code is causing a miscount or what
There is some work in progress here
https://github.com/bigdatagenomics/utils/pull/148/commits
https://github.com/bigdatagenomics/utils/pull/148/commits
ADAM version 0.31.0 will be released soon, with deprecations from utils-metrics 0.2.16 and others. Then ADAM version 0.32.0 will depend on utils version 0.3.0 and will also cross-build Spark 2 and Spark 3.
And I don't have access to any Spark 3 cluster environments for testing, e.g. AWS EMR is still at Spark 2.4.4
Block-gzipped (bgz/bgzf) FASTQ files are supported by the code in this package https://github.com/bigdatagenomics/adam/tree/master/adam-core/src/main/java/org/bdgenomics/adam/io
Most of it came in on this commit bigdatagenomics/adam@985e5d8
