<jupyter_start><jupyter_text>ChIP-Seq and RNA-SeqWe now want to combine our ChIP-Seq and RNA-Seq data, using pandas.<jupyter_code>import pandas as pd<jupyter_output><empty_output><jupyter_text>Gene annotationsWe can read our GFF file of gene annotations. It's just a text table.There are no headers, though. The columns of the file are:```['seqid', 'source', 'type', 'start', 'end', 'score', 'strand', 'phase', 'attributes']```<jupyter_code>genes = pd.read_csv("~/shared/MCB280A_data/S288C_R64-3-1/saccharomyces_cerevisiae_R64-3-1_20210421.gff",
                    delimiter='\t', header=None,
                   names=['seqid', 'source', 'type', 'start', 'end',
                          'score', 'strand', 'phase', 'attributes'])
genes.head()<jupyter_output><empty_output><jupyter_text>There are many different "types" of entries in the file. We can use the `value_counts()` method on the `type` column to see all of them.<jupyter_code>genes['type'].value_counts()<jupyter_output><empty_output><jupyter_text>We only want the `"gene"` type of entry and so we'll pick out just these rows out of the data frame.<jupyter_code>genes = genes[genes['type'] == "gene"]
genes.head()<jupyter_output><empty_output><jupyter_text>The systematic name of the gene is buried inside the `attributes` column. This column has a bunch of labeled data in the form _key_=_value_.We can extract it using the `.str.split()` method and then use `.str.replace()` to get rid of the wanted `ID=` part of it to create a `name` column.<jupyter_code>genes['name'] = genes['attributes'].str.split(';').str[0].str.replace("ID=", "")
genes.head()<jupyter_output><empty_output><jupyter_text>The data table lists the coordinates of genes, but we want to look for ChIP-Seq peaks in the promoter regions. For a `+` strand gene, the promoter is to the left (smaller coordinate numbers) and for a `-` strand gene it is to the right.We'll use a 1 kb window here.For a `+` strand gene, the promoter is _start_-1000 to _start_-1.For a `-` strand gene, the promoter is _end_+1 to _end_+1000.We can use the `np.where(...)` function from `numpy` to handle this situation. We'll compute the starting position for each promoter in the `prmstart` column.<jupyter_code>import numpy as np
genes['prmstart'] = np.where(genes['strand'] == "+", 
         genes['start'] - 1000,
         genes['end'] + 1)
genes.head()<jupyter_output><empty_output><jupyter_text>We'll be sure to check our coordinate math for genes on each strand.Next, we'll make a `prmend` column that is 1000 bases after the start. We don't need to do anything different based on the strand here.<jupyter_code>genes['prmend'] = genes['prmstart'] + 999
genes.head()<jupyter_output><empty_output><jupyter_text>ChIP-Seq PeaksNow, we'll read in the table of ChIP-Seq peaks.<jupyter_code>peaks = pd.read_csv("~/Hsf1/ChIP-Seq/macs2/Hsf1_ChIP_heatshk_peaks.xls",
                    delimiter = '\t', comment = '#')
peaks.head()<jupyter_output><empty_output><jupyter_text>At this point we need to loop over each peak and find out which promoter(s) it affects.We'll do this in a multi-step process.1. Generate a dictionary where keys are gene names and values are associated ChIP-Seq peak names. There won't be an entry for every gene.1. Convert this dictionary into a pandas `Series` and merge it in to a new column of the genes data frame, holding peak names. It will have many `NaN` entries for "missing" data.1. Merge the whole table of peaks (with enrichment and p-values) into the gene table based on the names.First, we will loop over each peak, using the `.itertuples()` method. In class, we'll just use the top 10 peaks by significance.<jupyter_code>gene_peaks = {}
top_peaks = peaks.sort_values('fold_enrichment', ascending=False).head(10)
for peak in top_peaks.itertuples():
    for gene in genes.itertuples():
        if (peak.chr == gene.seqid) and (peak.abs_summit >= gene.prmstart) and (peak.abs_summit <= gene.prmend):
            gene_peaks[gene.name] = peak.name
gene_peaks<jupyter_output><empty_output><jupyter_text>Now we need to convert this into a `Series` and give it a name.<jupyter_code>gene_peaks = pd.Series(gene_peaks, name='peak')
gene_peaks<jupyter_output><empty_output><jupyter_text>We can merge our named `Series` into the data frame of genes. We want to match up the `name` column in the genes table with the "index" of the `Series.<jupyter_code>genes2 = pd.merge(genes, gene_peaks,
                  left_on='name', right_index=True, how='left')
genes2.head()<jupyter_output><empty_output><jupyter_text>To make sure this worked, check on the row for the gene `YAL005C`, which does have a peak.<jupyter_code>genes2[genes2['name'] == "YGR210C"]<jupyter_output><empty_output><jupyter_text>Now, we will merge in the peaks table by matching up the `peak` column with the `name` column in the peaks table.<jupyter_code>genes3 = pd.merge(genes2, peaks,
                  left_on='peak', right_on='name', how='left')
genes3.head()<jupyter_output><empty_output><jupyter_text>Again, we'll need to check on a row that has a peak to be sure it worked.<jupyter_code>genes3[genes3['name_x'] == "YGR210C"]<jupyter_output><empty_output><jupyter_text>RNA-Seq dataFinally, we're ready to read in the table of RNA-Seq results.<jupyter_code>results = pd.read_csv("~/Hsf1/RNA-Seq/full.results.csv",
                      index_col = 0)
results.head()<jupyter_output><empty_output><jupyter_text>Again we can pull out some statistically significant genes to be sure it worked. Notice that almost every significant gene is down-regulated.<jupyter_code>results[results['padj'] < 0.1]<jupyter_output><empty_output><jupyter_text>Finally, we're ready to merge the results with the genes.We'll match up the `name_x` column of the genes (which was renamed because the peak table also had a name column) with the "index" of the results table.<jupyter_code>genes4 = pd.merge(genes3, results,
                  left_on='name_x', right_index=True, how='left')
genes4.head()<jupyter_output><empty_output><jupyter_text>Let's check out the ChIP-Seq genes in the RNA-Seq data.<jupyter_code>genes4.sort_values('fold_enrichment', ascending=False).head()<jupyter_output><empty_output><jupyter_text>Of course, most genes have a `NaN` missing value for fold enrichment. We can use the `pd.isna()` function to test whether a value is `NaN` or not.Rows where `fold_enrichment` is _not_ `NaN` are genes with a potential ChIP-Seq peak.<jupyter_code>pd.isna(genes4['fold_enrichment']).value_counts()<jupyter_output><empty_output><jupyter_text>Many of these genes have significant adjusted p-values and negative log fold-changes. Let's look at this trend more rigorously by plotting the histogram of fold-changes for these groups.Start by importing matplotlib.pyplot.<jupyter_code>import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>We'll make a histogram of `fold_enrichment` values for all genes.<jupyter_code>plt.hist(genes4['log2FoldChange'])<jupyter_output><empty_output><jupyter_text>Then, we'll make a similar histogram, for genes with a ChIP-Seq peak.We'll use the `range` parameter so the histograms are easier to compare.<jupyter_code>genes_with_peak = genes4[~pd.isna(genes4['fold_enrichment'])]
plt.hist(genes_with_peak['log2FoldChange'], range=(-4,4))<jupyter_output><empty_output><jupyter_text>This looks different—but is it significant?We'll use the `mannwhitneyu` function from `scipy.stats` to run a statistical test. To do this, we'll need to remove the `NaN` values using the `.dropna()` method.<jupyter_code>import scipy.stats as stats

nochip = genes4[pd.isna(genes4['fold_enrichment'])]['log2FoldChange'].dropna()
yeschip = genes4[~pd.isna(genes4['fold_enrichment'])]['log2FoldChange'].dropna()

stats.mannwhitneyu(nochip, yeschip)<jupyter_output><empty_output><jupyter_text>Finally, we'll make a table of the high-confidence targets that have a ChIP-Seq peak and an expression change.<jupyter_code>targets = genes4[(~pd.isna(genes4['fold_enrichment']))
                 & (genes4['padj'] < 0.1)]
targets<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11594202898550725, 'fraction_numerical': 0.012510838597795119, 'mean_word_length': 5.423229912490056, 'pattern_counts': {'":': 0, '<': 100, '<?xml version=': 0, '>': 98, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '24', 'dir': 'jupyter-structured-clean-dedup', 'id': '454257', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 2799, 'n_tokens_neox': 2592, 'n_words': 818, 'path': '09_multi_omics/09_a_genes_to_peaks from class.ipynb', 'repo_name': 'ingolia/mcb280-2021', 'size': '8073'}