Skip to content
Snippets Groups Projects
Commit 144bd3f4 authored by Beibei Chen's avatar Beibei Chen
Browse files

modify runMEME and set default for toppeak

parent b585c2ff
No related merge requests found
...@@ -114,7 +114,7 @@ workflow_parameters: ...@@ -114,7 +114,7 @@ workflow_parameters:
required: true required: true
description: | description: |
The number of top peaks to use for motif discovery. The number of top peaks to use for motif discovery.
default: 200
# ----------------------------------------------------------------------------- # -----------------------------------------------------------------------------
# SHINY APP CONFIGURATION # SHINY APP CONFIGURATION
......
...@@ -9,7 +9,6 @@ import string ...@@ -9,7 +9,6 @@ import string
import argparse as ap import argparse as ap
import logging import logging
import twobitreader import twobitreader
import pybedtools
import subprocess import subprocess
import pandas as pd import pandas as pd
from Bio import SeqIO from Bio import SeqIO
...@@ -52,7 +51,7 @@ def run_wrapper(args): ...@@ -52,7 +51,7 @@ def run_wrapper(args):
def run(infile, genome, limit, output): def run(infile, genome, limit, output):
infile = pybedtools.BedTool(infile) infile = open(infile).readlines()
logging.debug(len(infile)) logging.debug(len(infile))
genome = twobitreader.TwoBitFile(genome) genome = twobitreader.TwoBitFile(genome)
outfile = open(output+".fa","w") outfile = open(output+".fa","w")
...@@ -65,19 +64,20 @@ def run(infile, genome, limit, output): ...@@ -65,19 +64,20 @@ def run(infile, genome, limit, output):
rowcount += 1 rowcount += 1
#logging.debug(record) #logging.debug(record)
if rowcount <=limit: if rowcount <=limit:
#logging.debug(rowcount) seqbuf = record.rstrip().split("\t")
try: try:
#logging.debug(record.chrom) #logging.debug(record.chrom)
seq = genome[record.chrom][record.start:record.stop] seq = genome[seqbuf[0]][int(seqbuf[1]):int(seqbuf[2])]
except: except:
pass pass
else: else:
if record.strand == "-": if len(seqbuf)>=5:
seq = rc(seq) if seqbuf[5] == "-":
if len(record.fields)>=4: seq = rc(seq)
newfa_name = record.name if len(seqbuf)>=4:
newfa_name = seqbuf[3]
else: else:
newfa_name = "_".join(record.fields) newfa_name = "_".join(seqbuf)
newfa = SeqRecord(Seq(seq),newfa_name,description="") newfa = SeqRecord(Seq(seq),newfa_name,description="")
#logging.debug(seq) #logging.debug(seq)
SeqIO.write(newfa,outfile,"fasta") SeqIO.write(newfa,outfile,"fasta")
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment