from htsworkflow.pipelines.sequences import \
create_sequence_table, \
scan_for_sequences
+from htsworkflow.pipelines import qseq2fastq
+from htsworkflow.pipelines import srf2fastq
def main(cmdline=None):
parser = make_parser()
for a in args:
library_result_map.extend(read_library_result_map(a))
+ if opts.make_tree_from is not None:
+ make_tree_from(opts.make_tree_from, library_result_map)
+
if opts.daf is not None:
link_daf(opts.daf, library_result_map)
parser = OptionParser()
# commands
+ parser.add_option('--make-tree-from',
+ help="create directories & link data files",
+ default=None)
parser.add_option('--fastq', help="generate scripts for making fastq files",
default=False, action="store_true")
return parser
+def make_tree_from(source_path, library_result_map):
+ """Create a tree using data files from source path.
+ """
+ for lib_id, lib_path in library_result_map:
+ if not os.path.exists(lib_path):
+ logging.info("Making dir {0}".format(lib_path))
+ os.mkdir(lib_path)
+ source_lib_dir = os.path.join(source_path, lib_path)
+ if os.path.exists(source_lib_dir):
+ pass
+ for filename in os.listdir(source_lib_dir):
+ source_pathname = os.path.join(source_lib_dir, filename)
+ target_pathname = os.path.join(lib_path, filename)
+ if not os.path.exists(source_pathname):
+ raise IOError("{0} does not exist".format(source_pathname))
+ if not os.path.exists(target_pathname):
+ os.symlink(source_pathname, target_pathname)
+ logging.info(
+ 'LINK {0} to {1}'.format(source_pathname, target_pathname))
+
def build_fastqs(host, apidata, sequences_path, library_result_map,
force=False ):
"""
"""
qseq_condor_header = """
Universe=vanilla
-executable=/woldlab/rattus/lvol0/mus/home/diane/proj/solexa/gaworkflow/scripts/qseq2fastq
+executable=%(exe)s
error=log/qseq2fastq.err.$(process).log
output=log/qseq2fastq.out.$(process).log
log=log/qseq2fastq.log
-"""
+""" % {'exe': sys.executable }
qseq_condor_entries = []
srf_condor_header = """
Universe=vanilla
-executable=/woldlab/rattus/lvol0/mus/home/diane/proj/solexa/gaworkflow/scripts/srf2fastq
+executable=%(exe)s
output=log/srf_pair_fastq.out.$(process).log
error=log/srf_pair_fastq.err.$(process).log
log=log/srf_pair_fastq.log
environment="PYTHONPATH=/home/diane/lib/python2.6/site-packages:/home/diane/proj/solexa/gaworkflow PATH=/woldlab/rattus/lvol0/mus/home/diane/bin:/usr/bin:/bin"
-"""
+""" % {'exe': sys.executable }
srf_condor_entries = []
lib_db = find_archive_sequence_files(host,
apidata,
# write other lines
submission_files = os.listdir(result_dir)
fastqs = {}
+ fastq_attributes = {}
for f in submission_files:
attributes = view_map.find_attributes(f, lib_id)
if attributes is None:
continue
elif attributes.get("type", None) == 'fastq':
fastqs.setdefault(ext, set()).add(f)
+ fastq_attributes[ext] = attributes
else:
inifile.extend(
make_submission_section(line_counter,
)
inifile += ['']
line_counter += 1
+ # add in fastqs on a single line.
- # add in fastqs on a single line.
- for extension, fastq_set in fastqs.items():
+ for extension, fastq_files in fastqs.items():
inifile.extend(
make_submission_section(line_counter,
- fastq_set,
- attributes[extension])
+ fastq_files,
+ fastq_attributes[extension])
)
inifile += ['']
line_counter += 1
script = """Universe = vanilla
Executable = /bin/tar
-arguments = czvf ../%(archivename)s %(filelist)s
+arguments = czvhf ../%(archivename)s %(filelist)s
Error = compress.err.$(Process).log
Output = compress.out.$(Process).log
Log = /tmp/submission-compress-%(user)s.log
initialdir = %(initialdir)s
+environment="GZIP=-3"
+request_memory = 20
queue
"""
def condor_srf_to_fastq(srf_file, target_pathname, paired, flowcell=None,
mid=None, force=False):
- args = [ srf_file, ]
+ py = srf2fastq.__file__
+ args = [ py, srf_file, ]
if paired:
args.extend(['--left', target_pathname])
# this is ugly. I did it because I was pregenerating the target
def condor_qseq_to_fastq(qseq_file, target_pathname, flowcell=None, force=False):
- args = ['-i', qseq_file, '-o', target_pathname ]
+ py = qseq2fastq.__file__
+ args = [py, '-i', qseq_file, '-o', target_pathname ]
if flowcell is not None:
args.extend(['-f', flowcell])
script = """
self.apidata = apidata
self.lib_cache = {}
+ self.lib_paired = {}
# ma is "map algorithm"
ma = 'TH1014'
self.patterns = [
('*.bai', None),
- ('*.bam', self._guess_bam_view),
('*.splices.bam', 'Splices'),
+ ('*.bam', self._guess_bam_view),
+ ('junctions.bed', 'Junctions'),
('*.jnct', 'Junctions'),
+ ('*.unique.bigwig', None),
('*.plus.bigwig', 'PlusSignal'),
('*.minus.bigwig', 'MinusSignal'),
('*.bigwig', 'Signal'),
('*.condor', None),
('*.daf', None),
('*.ddf', None),
- ('cufflinks-0.9.0-genes.expr', 'GeneDeNovo'),
- ('cufflinks-0.9.0-transcripts.expr', 'TranscriptDeNovo'),
- ('cufflinks-0.9.0-transcripts.gtf', 'GeneModel'),
- ('GENCODE-v3c-genes.expr', 'GeneGencV3c'),
- ('GENCODE-v3c-transcripts.expr', 'TranscriptGencV3c'),
- ('GENCODE-v4-genes.expr', 'GeneGencV4'),
- ('GENCODE-v4-transcripts.expr', 'TranscriptGencV4'),
- ('GENCODE-v4-transcript.expr', 'TranscriptGencV4'),
+ ('*.?ufflinks-0.9.0?genes.expr', 'GeneDeNovo'),
+ ('*.?ufflinks-0.9.0?transcripts.expr', 'TranscriptDeNovo'),
+ ('*.?ufflinks-0.9.0?transcripts.gtf', 'GeneModel'),
+ ('*.GENCODE-v3c?genes.expr', 'GeneGCV3c'),
+ ('*.GENCODE-v3c?transcript*.expr', 'TranscriptGCV3c'),
+ ('*.GENCODE-v3c?transcript*.gtf', 'TranscriptGencV3c'),
+ ('*.GENCODE-v4?genes.expr', None), #'GeneGCV4'),
+ ('*.GENCODE-v4?transcript*.expr', None), #'TranscriptGCV4'),
+ ('*.GENCODE-v4?transcript*.gtf', None), #'TranscriptGencV4'),
+ ('*_1.75mers.fastq', 'FastqRd1'),
+ ('*_2.75mers.fastq', 'FastqRd2'),
('*_r1.fastq', 'FastqRd1'),
('*_r2.fastq', 'FastqRd2'),
('*.fastq', 'Fastq'),
('*.gtf', 'GeneModel'),
('*.ini', None),
('*.log', None),
- ('*.stats.txt', 'InsLength'),
+ ('paired-end-distribution*', 'InsLength'),
+ ('*.stats.txt', 'InsLength'),
('*.srf', None),
('*.wig', None),
('*.zip', None),
self.views = {
None: {"MapAlgorithm": "NA"},
"Paired": {"MapAlgorithm": ma},
+ "Aligns": {"MapAlgorithm": ma},
"Single": {"MapAlgorithm": ma},
"Splices": {"MapAlgorithm": ma},
"Junctions": {"MapAlgorithm": ma},
"GeneModel": {"MapAlgorithm": ma},
"GeneDeNovo": {"MapAlgorithm": ma},
"TranscriptDeNovo": {"MapAlgorithm": ma},
- "GeneGencV3c": {"MapAlgorithm": ma},
+ "GeneGCV3c": {"MapAlgorithm": ma},
+ "TranscriptGCV3c": {"MapAlgorithm": ma},
"TranscriptGencV3c": {"MapAlgorithm": ma},
- "GeneGencV4": {"MapAlgorithm": ma},
- "TranscriptGencV4": {"MapAlgorithm": ma},
+ "GeneGCV4": {"MapAlgorithm": ma},
+ "TranscriptGCV4": {"MapAlgorithm": ma},
"FastqRd1": {"MapAlgorithm": "NA", "type": "fastq"},
"FastqRd2": {"MapAlgorithm": "NA", "type": "fastq"},
"Fastq": {"MapAlgorithm": "NA", "type": "fastq" },
- "GeneModel": {"MapAlgorithm": ma},
"InsLength": {"MapAlgorithm": ma},
}
# view name is one of the attributes
'cell': lib_info['cell_line'],
'replicate': lib_info['replicate'],
}
- is_paired = self._is_paired(lib_info)
+ is_paired = self._is_paired(lib_id, lib_info)
if is_paired:
attributes.update(self.get_paired_attributes(lib_info))
if is_paired:
return "Paired"
else:
- return "Align"
+ return "Aligns"
- def _is_paired(self, lib_info):
+ def _is_paired(self, lib_id, lib_info):
"""Determine if a library is paired end"""
if len(lib_info["lane_set"]) == 0:
return False
-
- is_paired = 0
- isnot_paired = 0
- # check to see if all the flowcells are the same.
- # otherwise we might need to do something complicated
- for flowcell in lib_info["lane_set"]:
- if flowcell["paired_end"]:
- is_paired += 1
- else:
- isnot_paired += 1
- logging.debug("Library %s: %d were, %d were not paired" % \
- (lib_info["library_id"], is_paired, isnot_paired))
+ if not self.lib_paired.has_key(lib_id):
+ is_paired = 0
+ isnot_paired = 0
+ failed = 0
+ # check to see if all the flowcells are the same.
+ # otherwise we might need to do something complicated
+ for flowcell in lib_info["lane_set"]:
+ # yes there's also a status code, but this comparison
+ # is easier to read
+ if flowcell["status"].lower() == "failed":
+ # ignore failed flowcell
+ failed += 1
+ pass
+ elif flowcell["paired_end"]:
+ is_paired += 1
+ else:
+ isnot_paired += 1
+
+ logging.debug("Library %s: %d paired, %d single, %d failed" % \
+ (lib_info["library_id"], is_paired, isnot_paired, failed))
- if is_paired > isnot_paired:
- return True
- elif is_paired < isnot_paired:
- return False
- else:
- raise RuntimeError("Assumptions about paired vs not paired are wrong")
+ if is_paired > isnot_paired:
+ self.lib_paired[lib_id] = True
+ elif is_paired < isnot_paired:
+ self.lib_paired[lib_id] = False
+ else:
+ raise RuntimeError("Equal number of paired & unpaired lanes."\
+ "Can't guess library paired status")
+
+ return self.lib_paired[lib_id]
def get_paired_attributes(self, lib_info):
if lib_info['insert_size'] is None:
projects / htsworkflow.git / blobdiff