<jupyter_start><jupyter_text>Hsf1 binding motifWe have a table of Hsf1 ChIP-Seq peaks, including the name of the chromosome and the position of the "summit" within the peak.In earlier classes, we learned how to read in the yeast genome sequence.Today, we'll look up the genome sequence around each ChIP-Seq summit and write our own Fasta-format file that we can use as input to MEME, the motif search tool. We'll use pandas quite a bit here.<jupyter_code>import pandas as pd<jupyter_output><empty_output><jupyter_text>We saw last time that the "peaks" file from MACS contained all the information we wanted, with nicely named columns.We'll read in the peaks file `ChIP_1M_peaks.xls` as a pandas `DataFrame` and store it in the `peaks` variable<jupyter_code>peaks = pd.read_csv("~/Hsf1/ChIP-Seq/macs2/ChIP_1M_peaks.xls",
                    comment='#', delimiter='\t')
peaks.head()<jupyter_output><empty_output><jupyter_text>Use `sort_values()` to sort peaks according to their significance. Start with the most significant peak first, using the `ascending` parameter.<jupyter_code>peaks = peaks.sort_values("-log10(qvalue)", ascending=False)<jupyter_output><empty_output><jupyter_text>We can use the `head` function to take only the "best" peaks. We'll test our code on the best 5 peaks from the mini ChIP-Seq data set here in class today.<jupyter_code>peaks.head(5)<jupyter_output><empty_output><jupyter_text>We will extract the DNA sequence in a somewhat arbitrary window of ~100 base pairs, 50 on each side of the actual summit.To do this, define two new columns in the `peaks` data frame, one containing the start of the motif search region and the other containing the end of the region. We'll call these `motif_start` and `motif_end`<jupyter_code>peaks["motif_start"] = peaks["abs_summit"] - 50
peaks["motif_end"] = peaks["abs_summit"] + 50
peaks.head()<jupyter_output><empty_output><jupyter_text>We want to _loop_ over each row of the `DataFrame`.The `itertuples()` method gives us an iterator over the rows. Each row is a `namedtuple` with named fields taken from the column names in the DataFrame.To illustrate this, loop over each peak and print the name and the location (chromosome and start to end coordinates).<jupyter_code>for peak in peaks.head().itertuples():
    print('Peak {} on chromosome {} from {} to {}'
          .format(peak.name, peak.chr, peak.motif_start, peak.motif_end))<jupyter_output>Peak ChIP_1M_peak_2 on chromosome chrI from 141723 to 141823
Peak ChIP_1M_peak_83 on chromosome chrXII from 490163 to 490263
Peak ChIP_1M_peak_48 on chromosome chrVII from 914970 to 915070
Peak ChIP_1M_peak_37 on chromosome chrVI from 210322 to 210422
Peak ChIP_1M_peak_45 on chromosome chrVII from 772084 to 772184<jupyter_text>To extract the sequences from these regions of the genome, we need to refer back to the yeast genome sequence.Below, we import `SeqIO` from the `Bio` package.I've also provided the name of the fasta file of chromosome sequences.```chrom_file = '/home/jovyan/shared/MCB280A_data/S288C_R64-3-1/S288C_reference_sequence_R64-3-1_20210421.fsa'```<jupyter_code>from Bio import SeqIO

chrom_file = '/home/jovyan/shared/MCB280A_data/S288C_R64-3-1/S288C_reference_sequence_R64-3-1_20210421.fsa'<jupyter_output><empty_output><jupyter_text>As a starting point, we'll extract the sequence for one _manually defined_ genomic region.```motif_chr = 'chrI'motif_start = 141723motif_end = 141823```We'll loop through each chromosome, test whether its name (the `.name` field) matches our chromosome of interest, and then slice out the desired region from the sequence (in the `.seq` field).<jupyter_code>motif_chr = 'chrI'
motif_start = 141723
motif_end = 141823

chroms = SeqIO.parse(chrom_file, "fasta")
for chrom in chroms:
    if chrom.name == motif_chr:
        site_seq = chrom.seq[motif_start:motif_end]
        print(str(site_seq))<jupyter_output>AAAATCTGGAGAAAAAACGCAAATGACAGCTTCTAAACGTTCCGTGTGCTTTCTTTCTAGAATGTTCTGGAAAGTTTACAACAATCCACAAGAACGAAAA<jupyter_text>Next, instead of manually defining the region, we'll loop over the rows in `peaks` and use the values in the table for each sequence.<jupyter_code>for peak in peaks.head().itertuples():
    chroms = SeqIO.parse(chrom_file, "fasta")
    for chrom in chroms:
        if chrom.name == peak.chr:
            site = chrom.seq[peak.motif_start:peak.motif_end]
            print('{}\t{}'.format(peak.name, site))<jupyter_output>ChIP_1M_peak_2	AAAATCTGGAGAAAAAACGCAAATGACAGCTTCTAAACGTTCCGTGTGCTTTCTTTCTAGAATGTTCTGGAAAGTTTACAACAATCCACAAGAACGAAAA
ChIP_1M_peak_83	AGGGTAAAACTAACCTGTCTCACGACGGTCTAAACCCAGCTCACGTTCCCTATTAGTGGGTGAACAATCCAACGCTTACCGAATTCTGCTTCGGTATGAT
ChIP_1M_peak_48	AAAAGCCCTTACTAACCCTACAATAAATTGTGCCGAAACCCTCTGGAGTTTTCTAGAATATTCTAGCCCCATCAGGGCTAGAATATCCTAAAAGTTTATA
ChIP_1M_peak_37	TCACGTTATCAGGCTCATAGCTTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGATTGTTGTTCTAGTCG
ChIP_1M_peak_45	GTTGTAAAGTTACTGACACTTTTTTTTCTAGAAAGTTCCGGAAAATTGCGACACTCGGTGGAGCTCGAGAGTTGTATCCAGTTTTCTTGTTCGGCGATAT<jupyter_text>This works okay, but it's inefficient because we read the genome file once for each peak. Instead, we can read the genome just once and store the sequence in a `dict`. We'll use the chromosome names as the _key_ and the sequence as the _value_.Build this dictionary in a variable named `genome`.<jupyter_code>genome = {}
chroms = SeqIO.parse(chrom_file, 'fasta')
for chrom in chroms:
    genome[chrom.name] = chrom.seq<jupyter_output><empty_output><jupyter_text>Now, we can use the `.get()` method on the `genome` dictionary to look up any chromosome sequence we want, very quickly.Adapt the loop above for this more efficient<jupyter_code>for peak in peaks.head().itertuples():
    chrom_seq = genome.get(peak.chr)
    site = chrom_seq[peak.motif_start:peak.motif_end]
    print('{}\t{}'.format(peak.name, site))

#for peak in peaks.head().itertuples():
#    chroms = SeqIO.parse(chrom_file, "fasta")
#    for chrom in chroms:
#        if chrom.name == peak.chr:
#            site = chrom.seq[peak.motif_start:peak.motif_end]
#            print('{}\t{}'.format(peak.name, site))<jupyter_output>ChIP_1M_peak_2	AAAATCTGGAGAAAAAACGCAAATGACAGCTTCTAAACGTTCCGTGTGCTTTCTTTCTAGAATGTTCTGGAAAGTTTACAACAATCCACAAGAACGAAAA
ChIP_1M_peak_83	AGGGTAAAACTAACCTGTCTCACGACGGTCTAAACCCAGCTCACGTTCCCTATTAGTGGGTGAACAATCCAACGCTTACCGAATTCTGCTTCGGTATGAT
ChIP_1M_peak_48	AAAAGCCCTTACTAACCCTACAATAAATTGTGCCGAAACCCTCTGGAGTTTTCTAGAATATTCTAGCCCCATCAGGGCTAGAATATCCTAAAAGTTTATA
ChIP_1M_peak_37	TCACGTTATCAGGCTCATAGCTTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGATTGTTGTTCTAGTCG
ChIP_1M_peak_45	GTTGTAAAGTTACTGACACTTTTTTTTCTAGAAAGTTCCGGAAAATTGCGACACTCGGTGGAGCTCGAGAGTTGTATCCAGTTTTCTTGTTCGGCGATAT<jupyter_text>Now, let's write these sequences into our own Fasta file.We'll need to create our own `Seq` and `SeqRecord` objects. To do this, we'll need to import some more parts of biopython:```from Bio.Seq import Seqfrom Bio.SeqRecord import SeqRecord```<jupyter_code>from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord<jupyter_output><empty_output><jupyter_text>After doing this, we can create our own sequence object by writing `Seq(...)`<jupyter_code>Seq("ACGTACGT")<jupyter_output><empty_output><jupyter_text>We can create our own sequence record—a sequence object, plus a _name_ in the `id` field, as well.<jupyter_code>SeqRecord(Seq("ACGTACGT"), id="MySeq", description='')<jupyter_output><empty_output><jupyter_text>Finally, we can use the `SeqIO.write()` function to write a list of `SeqRecord`s into a file. This function has three parameters:```SeqIO.write(list-of-sequences, filename, format)```<jupyter_code>rec1 = SeqRecord(Seq("ACGTACGT"), id="MySeq", description='')
rec2 = SeqRecord(Seq("GCATGCAT"), id="AnotherSeq", description='another sequence')
SeqIO.write([rec1, rec2], "test.fa", "fasta")<jupyter_output><empty_output><jupyter_text>Now, we will construct a list called `sites` containing a `SeqRecord` for each ChIP-Seq peak. Then, we'll write the list into a file called `"Hsf1_sites.fa"`.<jupyter_code>site_seqs = []

for peak in peaks.head().itertuples():
    chrom_seq = genome.get(peak.chr)
    site = chrom_seq[peak.motif_start:peak.motif_end]
    description = '{}:{}-{}'.format(peak.chr, 
                                    peak.motif_start, 
                                    peak.motif_end)
    rec = SeqRecord(site, id=peak.name, description=description)
    site_seqs.append(rec)

# STUFF HERE

SeqIO.write(site_seqs, "Hsf1_sites.fa", "fasta")
peaks.head(10)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08768022506154026, 'fraction_numerical': 0.023795569100926034, 'mean_word_length': 6.022222222222222, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 61, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '642224', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 3320, 'n_tokens_neox': 2957, 'n_words': 772, 'path': '05_chip_seq_motifs/01_extract_seqs from class.ipynb', 'repo_name': 'ingolia/mcb280-2021', 'size': '8531'}