<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Read motifs from files in other formats.
Step1:  You can convert a motif to several formats.
Step2:  Some other useful tidbits.
Step3:  To convert a motif to an image, use to_img(). Supported formats are png, ps and pdf.
Step4:  Motif scanning
For very simple scanning, you can just use a Motif instance. Let’s say we have a FASTA file called test.small.fa that looks like this
Step5:  This return a dictionary with the sequence names as keys. The value is a list with positions where the motif matches. Here, as the AP1 motif is a palindrome, you see matches on both forward and reverse strand. This is more clear when we use pwm_scan_all() that returns position, score and strand for every match.
Step6:  The number of matches to return is set to 50 by default, you can control this by setting the nreport argument. Use scan_rc=False to only scan the forward orientation.
Step7:  While this functionality works, it is not very efficient. To scan many motifs in potentially many sequences, use the functionality in the scanner module. If you only want the best match per sequence, there is a utility function called scan_to_best_match, otherwise, use the Scanner class.
Step8:  The matches are in the same order as the sequences in the original file.
While this function can be very useful, a Scanner instance is much more flexible. You can scan different input formats (BED, FASTA, regions), and control the thresholds and output.
As an example we will use the file Gm12878.CTCF.top500.w200.fa that contains 500 top CTCF peaks. We will get the CTCF motif and scan this file in a number of different ways.
Step9:  Now let’s get the best score for the CTCF motif for each sequence.
Step10:  In many cases you’ll want to set a threshold. In this example we’ll use a 1% FPR threshold, based on scanning randomly selected sequences from the hg38 genome. The first time you run this, it will take a while. However, the tresholds will be cached. This means that for the same combination of motifs and genome, the previously generated threshold will be used.
Step11:  Finding de novo motifs
Let’s take the Gm12878.CTCF.top500.w200.fa file as example again. For a basic example we’ll just use two motif finders, as they’re quick to run.
Step12:  This will basically run the same pipeline as the gimme motifs command. All output files will be stored in outdir and gimme_motifs returns a list of Motif instances. If you only need the motifs but not the graphical report, you can decide to skip it by setting create_report to False. Additionally, you can choose to skip clustering (cluster=False) or to skip calculation of significance (filter_significant=False). For instance, the following command will only predict motifs and cluster them.
Step13:  All parameters for motif finding are set by the params argument
Although the gimme_motifs function is probably the easiest way to run the de novo finding tools, you can also run any of the tools directly. In this case you would also have to supply the background file if the specific tool requires it.
Step14:  Motif statistics
With some motifs, a sample file and a background file you can calculate motif statistics. Let’s say I wanted to know which of the p53-family motifs is most enriched in the file TAp73alpha.fa.
First, we’ll generate a GC%-matched genomic background. Then we only select p53 motifs.
Step15:  A lot of statistics are generated and you will not always need all of them. You can choose one or more specific metrics with the additional stats argument.
Step16:  Motif comparison
Step17:  Compare two motifs
Step18:  Find closest match in a motif database
<ASSISTANT_TASK:>
Python Code:
with open("MA0099.3.jaspar") as f:
    motifs = read_motifs(f, fmt="jaspar")
print(motifs[0])
Explanation: Read motifs from files in other formats.
End of explanation
with open("example.pfm") as f:
    motifs = read_motifs(f)
# pwm
print(motifs[0].to_pwm())
# pfm
print(motifs[0].to_pfm())
# consensus sequence
print(motifs[0].to_consensus())
# TRANSFAC
print(motifs[0].to_transfac())
# MEME
print(motifs[0].to_meme())
Explanation: You can convert a motif to several formats.
End of explanation
m = motif_from_consensus("NTGASTCAN")
print(len(m))
# Trim by information content
m.trim(0.5)
print(m.to_consensus(), len(m))
# Slices
print(m[:3].to_consensus())
# Shuffle
random_motif = motif_from_consensus("NTGASTGAN").randomize()
print(random_motif)
Explanation: Some other useful tidbits.
End of explanation
m = motif_from_consensus("NTGASTCAN")
m.to_img("ap1.png", fmt="png")
from IPython.display import Image
Image("ap1.png")
Explanation: To convert a motif to an image, use to_img(). Supported formats are png, ps and pdf.
End of explanation
from gimmemotifs.motif import motif_from_consensus
from gimmemotifs.fasta import Fasta
f = Fasta("test.small.fa")
m = motif_from_consensus("TGAsTCA")
m.pwm_scan(f)
Explanation: Motif scanning
For very simple scanning, you can just use a Motif instance. Let’s say we have a FASTA file called test.small.fa that looks like this:
```
seq1
AAAAAAAAAAAAAAAAAAAAAA
seq2
CGCGCGTGAGTCACGCGCGCGCG
seq3
TGASTCAAAAAAAAAATGASTCA
```
Now we can use this file for scanning.
End of explanation
m.pwm_scan_all(f)
Explanation: This return a dictionary with the sequence names as keys. The value is a list with positions where the motif matches. Here, as the AP1 motif is a palindrome, you see matches on both forward and reverse strand. This is more clear when we use pwm_scan_all() that returns position, score and strand for every match.
End of explanation
m.pwm_scan_all(f, nreport=1, scan_rc=False)
Explanation: The number of matches to return is set to 50 by default, you can control this by setting the nreport argument. Use scan_rc=False to only scan the forward orientation.
End of explanation
from gimmemotifs.motif import motif_from_consensus
from gimmemotifs.scanner import scan_to_best_match
m1 = motif_from_consensus("TGAsTCA")
m1.id = "AP1"
m2 = motif_from_consensus("CGCG")
m2.id = "CG"
motifs = [m1, m2]
print("motif\tpos\tscore")
result = scan_to_best_match("test.small.fa", motifs)
for motif, matches in result.items():
    for match in matches:
        print("{}\t{}\t{}".format(motif, match[1], match[0]))
Explanation: While this functionality works, it is not very efficient. To scan many motifs in potentially many sequences, use the functionality in the scanner module. If you only want the best match per sequence, there is a utility function called scan_to_best_match, otherwise, use the Scanner class.
End of explanation
from gimmemotifs.motif import default_motifs
from gimmemotifs.scanner import Scanner
from gimmemotifs.fasta import Fasta
import numpy as np
# Input file
fname = "Gm12878.CTCF.top500.w200.fa"
# Select the CTCF motif from the default motif database
motifs = [m for m in default_motifs() if "CTCF" in m.factors['direct']][:1]
# Initialize the scanner
s = Scanner()
s.set_motifs(motifs)
Explanation: The matches are in the same order as the sequences in the original file.
While this function can be very useful, a Scanner instance is much more flexible. You can scan different input formats (BED, FASTA, regions), and control the thresholds and output.
As an example we will use the file Gm12878.CTCF.top500.w200.fa that contains 500 top CTCF peaks. We will get the CTCF motif and scan this file in a number of different ways.
End of explanation
scores = [r[0] for r in s.best_score("Gm12878.CTCF.top500.w200.fa")]
print("{}\t{:.2f}\t{:.2f}\t{:.2f}".format(
    len(scores),
    np.mean(scores),
    np.min(scores),
    np.max(scores)
    ))
Explanation: Now let’s get the best score for the CTCF motif for each sequence.
End of explanation
# Set a 1% FPR threshold based on random hg38 sequence
s.set_genome("hg38")
s.set_threshold(fpr=0.01)
# get the number of sequences with at least one match
counts = [n[0] for n in s.count("Gm12878.CTCF.top500.w200.fa", nreport=1)]
print(counts[:10])
# or the grand total of number of sequences with 1 match
print(s.total_count("Gm12878.CTCF.top500.w200.fa", nreport=1))
# Scanner.scan() just gives all information
seqs = Fasta("Gm12878.CTCF.top500.w200.fa")[:10]
for i,result in enumerate(s.scan(seqs)):
    seqname = seqs.ids[i]
    for m,matches in enumerate(result):
        motif = motifs[m]
        for score, pos, strand in matches:
            print(seqname, motif, score, pos, strand)
Explanation: In many cases you’ll want to set a threshold. In this example we’ll use a 1% FPR threshold, based on scanning randomly selected sequences from the hg38 genome. The first time you run this, it will take a while. However, the tresholds will be cached. This means that for the same combination of motifs and genome, the previously generated threshold will be used.
End of explanation
from gimmemotifs.denovo import gimme_motifs
peaks = "Gm12878.CTCF.top500.w200.fa"
outdir = "CTCF.gimme"
params = {
    "tools": "Homer,BioProspector",
    "genome": "hg38",
    }
motifs = gimme_motifs(peaks, outdir, params=params)
Explanation: Finding de novo motifs
Let’s take the Gm12878.CTCF.top500.w200.fa file as example again. For a basic example we’ll just use two motif finders, as they’re quick to run.
End of explanation
motifs = gimme_motifs(peaks, outdir,
    params=params, filter_significant=False, create_report=False)
Explanation: This will basically run the same pipeline as the gimme motifs command. All output files will be stored in outdir and gimme_motifs returns a list of Motif instances. If you only need the motifs but not the graphical report, you can decide to skip it by setting create_report to False. Additionally, you can choose to skip clustering (cluster=False) or to skip calculation of significance (filter_significant=False). For instance, the following command will only predict motifs and cluster them.
End of explanation
from gimmemotifs.tools import get_tool
from gimmemotifs.background import MatchedGcFasta
m = get_tool("homer")  # tool name is case-insensitive
# Create a background fasta file with a similar GC%
fa = MatchedGcFasta("TAp73alpha.fa", number=1000)
fa.writefasta("bg.fa")
# Run motif prediction
params = {
    "background": "bg.fa",
    "width": "20",
    "number": 5,
}
motifs, stdout, stderr = m.run("TAp73alpha.fa", params=params)
print(motifs[0].to_consensus())
Explanation: All parameters for motif finding are set by the params argument
Although the gimme_motifs function is probably the easiest way to run the de novo finding tools, you can also run any of the tools directly. In this case you would also have to supply the background file if the specific tool requires it.
End of explanation
from gimmemotifs.background import MatchedGcFasta
from gimmemotifs.fasta import Fasta
from gimmemotifs.stats import calc_stats
from gimmemotifs.motif import default_motifs
sample = "TAp73alpha.fa"
bg = MatchedGcFasta(sample, genome="hg19", number=1000)
motifs = [m for m in default_motifs() if any(f in m.factors['direct'] for f in ["TP53", "TP63", "TP73"])]
stats = calc_stats(motifs, sample, bg)
print("Stats for", motifs[0])
for k, v in stats[str(motifs[0])].items():
    print(k,v)
print()
best_motif = sorted(motifs, key=lambda x: stats[str(x)]["recall_at_fdr"])[-1]
print("Best motif (recall at 10% FDR):", best_motif)
Explanation: Motif statistics
With some motifs, a sample file and a background file you can calculate motif statistics. Let’s say I wanted to know which of the p53-family motifs is most enriched in the file TAp73alpha.fa.
First, we’ll generate a GC%-matched genomic background. Then we only select p53 motifs.
End of explanation
metrics = ["roc_auc", "recall_at_fdr"]
stats = calc_stats(motifs, sample, bg, stats=metrics)
for metric in metrics:
    for motif in motifs:
        print("{}\t{}\t{:.2f}".format(
            motif.id, metric, stats[str(motif)][metric]
            ))
Explanation: A lot of statistics are generated and you will not always need all of them. You can choose one or more specific metrics with the additional stats argument.
End of explanation
from gimmemotifs.comparison import MotifComparer
from gimmemotifs.motif import motif_from_consensus
from gimmemotifs.motif import read_motifs
Explanation: Motif comparison
End of explanation
m1 = motif_from_consensus("RRRCATGYYY")
m2 = motif_from_consensus("TCRTGT")
mc = MotifComparer()
score, pos, orient = mc.compare_motifs(m1, m2)
if orient == -1:
    m2 = m2.rc()
pad1, pad2 = "", ""
if pos < 0:
    pad1 = " " * -pos 
elif pos > 0:
    pad2 =" " * pos
print(pad1 + m1.to_consensus())
print(pad2 + m2.to_consensus())
Explanation: Compare two motifs
End of explanation
motifs = [
    motif_from_consensus("GATA"),
    motif_from_consensus("NTATAWA"),
    motif_from_consensus("ACGCG"),
]
mc = MotifComparer()
results = mc.get_closest_match(motifs, dbmotifs=read_motifs("HOMER"), metric="seqcor")
# Load motifs
db = read_motifs("HOMER", as_dict=True)
for motif in motifs:
    match, scores = results[motif.id]
    print("{}: {} - {:.3f}".format(motif.id, match, scores[0]))
    dbmotif = db[match]
    orient = scores[2]
    if orient == -1:
        dbmotif = dbmotif.rc()
    padm, padd = 0, 0
    if scores[1] < 0:
        padm = -scores[1]
    elif scores[1] > 0:
        padd = scores[1]
    print(" " * padm + motif.to_consensus())
    print(" " * padd + dbmotif.to_consensus())
    print()
Explanation: Find closest match in a motif database
End of explanation
<END_TASK>