<jupyter_start><jupyter_text>Exploratory Data Analysis of Genomic Datasets with ADAM and Mango  Configuring ADAM and Mango on EMR Mango uses docker containers to be run easily on EMR. To get everything setup and installed, follow EMR documentation at http://bdg-mango.readthedocs.io/en/latest/cloud/emr.html.  Loading Data from the 1000 Genomes ProjectIn this tutorial, we will use ADAM and Mango to discover interesting variants in the child of a 1000 Genomes trio.First, let’s import ADAM and Mango modules, as well as any Spark modules we need:<jupyter_code># Import ADAM modules
from bdgenomics.adam.adamContext import ADAMContext
from bdgenomics.adam.rdd import AlignmentDataset, CoverageDataset
from bdgenomics.adam.stringency import LENIENT, _toJava

# Import Mango modules
from bdgenomics.mango.alignments import *
from bdgenomics.mango.coverage import CoverageDistribution

# Import Spark modules
from pyspark.sql import functions as sf<jupyter_output><empty_output><jupyter_text>Next, we will create an ADAMContext. ADAMContext allows us to load and manipulate genomic data.<jupyter_code># Create ADAM Context
ac = ADAMContext(spark)<jupyter_output><empty_output><jupyter_text>Variant Analysis with Spark SQL In this analysis, we will view a trio (NA19685, NA19661, and NA19660) and search for variants that are present in the child but not present in the parents. These are interesting regions, as they may indicate sights of de novo variation that may contribute to multiple disorders.First, we will load in a subset of variant data from chromosome 17:<jupyter_code>pathPrefix = 's3://1000genomes/phase1/analysis_results/integrated_call_sets/'

genotypesPath = pathPrefix + 'ALL.chr17.integrated_phase1_v3.20101123.snps_indels_svs.genotypes.vcf.gz'
genotypes = ac.loadGenotypes(genotypesPath)

genotypes_df  = genotypes.toDF()<jupyter_output><empty_output><jupyter_text>We can take a look at the schema by printing the columns in the dataframe.<jupyter_code># cache genotypes and show the schema
genotypes_df.columns<jupyter_output><empty_output><jupyter_text>This genotypes dataset contains all samples from the 1000 Genomes Project. Therefore, we will next filter genotypes to only consider samples that are in the NA19685 trio.<jupyter_code># trio IDs
IDs = ['NA19685','NA19661','NA19660']

trio_df = genotypes_df.filter(genotypes_df["sampleId"].isin(IDs))<jupyter_output><empty_output><jupyter_text>We will next add a new column to our dataframe that determines the genomic location of each variant. This is defined by the chromosome (referenceName) and the start and end position of the variant.<jupyter_code># Add ReferenceRegion column and group by referenceRegion
trios_with_referenceRegion = trio_df.withColumn('ReferenceRegion', 
                    sf.concat(sf.col('referenceName'),sf.lit(':'), sf.col('start'), sf.lit('-'), sf.col('end')))<jupyter_output><empty_output><jupyter_text>Now, we want to query our dataset to find de novo variants. But first, we must register our dataframe with Spark SQL.<jupyter_code>#  Register df with Spark SQL
trios_with_referenceRegion.createOrReplaceTempView("trios")<jupyter_output><empty_output><jupyter_text>Now that our dataframe is registered, we can run SQL queries on it. For our first query, we will select the names of a subset of variants belonging to sample NA19685 that have at least one alternative (ALT) allele.<jupyter_code># filter by alleles. This is a list of variant names that have an alternate allele for the child
alternate_variant_sites = spark.sql("SELECT variant.names[0] AS snp FROM trios \
                                    WHERE array_contains(alleles, 'ALT') AND sampleId == 'NA19685'") 

collected_sites = list(map(lambda x: x.snp, alternate_variant_sites.take(100)))<jupyter_output><empty_output><jupyter_text>For our next query, we will filter a subset of sites in which the parents have both reference alleles. We then filter these variants by the set produced above from the child.<jupyter_code># get parent records and filter by only REF locations for variant names that were found in the child with an ALT
filtered1 = spark.sql("SELECT * FROM trios WHERE sampleId == 'NA19661' or sampleId == 'NA19660' \
            AND !array_contains(alleles, 'ALT')")


filtered2 = filtered1.filter(filtered1["variant.names"][0].isin(collected_sites))

snp_counts = filtered2.take(100)
# view snp names as a list
set([x.variant.names[0] for x in snp_counts])<jupyter_output><empty_output><jupyter_text>Working with Alignment Data Now, we can explore these specific variant sites in the raw genomic alignment data. First, let’s load in the data for the NA19685 trio:<jupyter_code># load in NA19685 exome from s3a
childReadsPath = 's3a://1000genomes/phase1/data/NA19685/exome_alignment/NA19685.mapped.illumina.mosaik.MXL.exome.20110411.bam'
parent1ReadsPath = 's3a://1000genomes/phase1/data/NA19660/exome_alignment/NA19660.mapped.illumina.mosaik.MXL.exome.20110411.bam'
parent2ReadsPath = 's3a://1000genomes/phase1/data/NA19661/exome_alignment/NA19661.mapped.illumina.mosaik.MXL.exome.20110411.bam'

childReads = ac.loadAlignments(childReadsPath, stringency=LENIENT)
parent1Reads = ac.loadAlignments(parent1ReadsPath, stringency=LENIENT)
parent2Reads = ac.loadAlignments(parent2ReadsPath, stringency=LENIENT)<jupyter_output><empty_output><jupyter_text>Quality Control of Alignment Data One popular analysis to visually re-affirm the quality of genomic alignment data is by viewing coverage distribution. Coverage distribution gives us an idea of the read coverage we have across a sample. Next, we will generate a sample coverage distribution plot for the child alignment data on chromosome 17.<jupyter_code># calculate read coverage
# Takes 2-3 minutes
childCoverage = childReads.transform(lambda x: x.filter(x.referenceName == "17")).toCoverage()<jupyter_output><empty_output><jupyter_text>Now that coverage data is calculated and cached, we will compute the coverage distribution of all three samples and plot the coverage distribution.<jupyter_code># Calculate coverage distribution
# You can check the progress in the SparkUI by navigating to 
# <PUBLIC_MASTER_DNS>:8088 and clicking on the currently running Spark application.
cd = CoverageDistribution(spark, childCoverage, bin_size = 1)
ax, results = cd.plotDistributions(normalize=True, cumulative=False)

ax.set_title("Coverage Distribution")
ax.set_ylabel("Counts")
ax.set_xlabel("Coverage Depth")
ax.set_xscale("log")
plt.show()<jupyter_output><empty_output><jupyter_text>Now that we are done with coverage, we can unpersist these datasets to clear space in memory for the next analysis.<jupyter_code>childCoverage.unpersist()<jupyter_output><empty_output><jupyter_text>Viewing Sites with Missense Variants in the ProbandAfter verifying alignment data and filtering variants, we have 4 genes with potential missense mutations in the proband, including YBX2, ZNF286B, KSR1, and GNA13. We can visually verify these sites by filtering and viewing the raw reads of the child and parents.First, let's view the child reads. If we zoom in to the location of the GNA13 variant (63052580-63052581) we can see a heterozygous T to A call.<jupyter_code># view missense variant at GNA13: 63052580-63052581 (SNP rs201316886) in child
# Takes about 2 minutes to collect data from workers
childViz = AlignmentSummary(spark, ac, childReads)
contig = "17"
start = 63052180
end = 63052981

childViz.viewPileup(contig, start, end)<jupyter_output><empty_output><jupyter_text>It looks like there indeed is a variant at this position, possibly a heterozygous SNP with alternate allele A. Let's look at the parent data to verify this variant does not appear in the parents.<jupyter_code># view missense variant at GNA13: 63052580-63052581 in parent 1
parent1Viz = AlignmentSummary(spark, ac, parent1Reads)
contig = "17"
start = 63052180
end = 63052981

parent1Viz.viewPileup(contig, start, end)
# view missense variant at GNA13: 63052580-63052581 in parent 2 
parent2Viz = AlignmentSummary(spark, ac, parent2Reads)
contig = "17"
start = 63052180
end = 63052981

parent2Viz.viewPileup(contig, start, end)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07468727005150846, 'fraction_numerical': 0.04329163600686779, 'mean_word_length': 6.024117140396211, 'pattern_counts': {'":': 0, '<': 62, '<?xml version=': 0, '>': 62, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '364261', 'license': 'ECL-2.0', 'n_tokens_mistral': 2787, 'n_tokens_neox': 2446, 'n_words': 861, 'path': 'example-files/notebooks/aws-1000genomes.ipynb', 'repo_name': 'heuermh/mango', 'size': '8154'}