Skip to content
Snippets Groups Projects
Commit 3d77cb5b authored by Venkat Malladi's avatar Venkat Malladi
Browse files

update the naming of the reads.

parent d14969c6
1 merge request!5Resolve "Add mapping and trimming"
Pipeline #1044 passed with stage
in 1 hour, 12 minutes, and 11 seconds
......@@ -112,24 +112,27 @@ def align_se(fastq, sai, reference, fastq_basename):
'''Use BWA to align SE data.'''
sam_filename = "%s.sam" % (fastq_basename)
bam_filename = '%s.bam' % (fastq_basename)
steps = [
"bwa samse %s %s %s"
% (reference, fastq[0], sai[0]),
"samtools view -@%d -Su -" % (cpu_count()),
"samtools sort -@%d -o %s"
% (cpu_count(), fastq_basename)]
% (cpu_count(), bam_filename)]
out, err = utils.run_pipe(steps)
if err:
logger.error("samse/samtools error: %s" % (err))
return bam_filename
def align_pe(fastq, sai, reference, fastq_basename):
'''Use BWA to align PE data.'''
sam_filename = "%s.sam" % (fastq_basename)
badcigar_filename = "%s.badReads" % (fastq_basename)
bam_filename = '%s.bam' % (fastq_basename)
# Remove read pairs with bad CIGAR strings and sort by position
steps = [
......@@ -150,12 +153,14 @@ def align_pe(fastq, sai, reference, fastq_basename):
"grep -v -F -f %s" % (badcigar_filename),
"samtools view -@%d -Su -" % (cpu_count()),
"samtools sort -@%d -o %s"
% (cpu_count(), fastq_basename)]
% (cpu_count(), bam_filename)]
out, err = utils.run_pipe(steps)
if err:
logger.error("samtools error: %s" % (err))
return bam_filename
def main():
args = get_args()
......@@ -184,18 +189,18 @@ def main():
strip_extensions(fastq[1], STRIP_EXTENSIONS))
fastq_basename = fastq_r1_basename + fastq_r2_basename
align_pe(fastq, sai, reference, fastq_basename)
bam_filename = align_pe(fastq, sai, reference, fastq_basename)
else:
fastq_basename = os.path.basename(
strip_extensions(fastq[0], STRIP_EXTENSIONS))
align_se(fastq, sai, reference, fastq_basename)
bam_filename = align_se(fastq, sai, reference, fastq_basename)
bam_mapstats_filename = '%s.raw.srt.bam.flagstat.qc' % (fastq_basename)
with open(bam_mapstats_filename, 'w') as fh:
subprocess.check_call(
shlex.split("samtools flagstat %s" % (bam_mapstats_filename)),
shlex.split("samtools flagstat %s" % (bam_filename)),
stdout=fh)
# Remove sai files
......
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