diff ICGC_STAR_ALIGNMENT_PIPELINE/star_align.py @ 0:0446d6f7d593 draft default tip

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ICGC_STAR_ALIGNMENT_PIPELINE/star_align.py	Wed Apr 18 02:38:42 2018 -0400
@@ -0,0 +1,527 @@
+#!/usr/bin/env python
+
+import os
+import sys
+import re
+import string
+import tempfile
+import subprocess
+import argparse
+import shutil
+import lxml.etree as etree
+import fnmatch
+
+def walk_dir(base, pattern):
+    files = []
+    for root, dirnames, filenames in os.walk(base):
+        for fname in fnmatch.filter(filenames, pattern):
+            files.append(os.path.join(root, fname))
+    return files
+
+def scan_workdir(base): 
+
+    ### scan for paired-end files
+    #############################
+
+    ### unzipped fastq input
+    fastq_files = walk_dir(base, "*_read[12]_*fastq")
+    if len(fastq_files):
+        o = {}
+        for i in sorted(fastq_files):
+            basename = re.sub(r'_read[12]', '', i)
+            try:
+                o[basename].append(i)
+            except KeyError:
+                o[basename] = [i]
+        if not all( (len(i) == 2 for i in o.values())):
+            raise Exception("Missing Pair")
+        return ( 'cat', list( (os.path.basename(i), o[i][0], o[i][1]) for i in o.keys()), 'PE') 
+
+    ### unzipped fastq input
+    fastq_files = walk_dir(base, "*_R[12]_001.fastq")
+    if len(fastq_files):
+        o = {}
+        for i in fastq_files:
+            basename = re.sub(r'_R[12]_001.fastq$', '', i)
+            o[basename] = o.get(basename, 0) + 1
+        if not all( (i == 2 for i in o.values())):
+            raise Exception("Missing Pair")
+        return ( 'cat', list( (os.path.basename(i), "%s_R1_001.fastq" % i,"%s_R2_001.fastq" % i) for i in o.keys()), 'PE') 
+
+    ### unzipped fastq input
+    fastq_files = walk_dir(base, "*_[12].fastq")
+    if len(fastq_files):
+        o = {}
+        for i in fastq_files:
+            basename = re.sub(r'_[12].fastq$', '', i)
+            o[basename] = o.get(basename, 0) + 1
+        if not all( (i == 2 for i in o.values())):
+            raise Exception("Missing Pair")
+        return ( 'cat', list( (os.path.basename(i), "%s_1.fastq" % i,"%s_2.fastq" % i) for i in o.keys()), 'PE') 
+
+    ### unzipped fastq input
+    fastq_files = walk_dir(base, "*[.][12].fastq")
+    if len(fastq_files):
+        o = {}
+        for i in fastq_files:
+            basename = re.sub(r'[.][12].fastq$', '', i)
+            o[basename] = o.get(basename, 0) + 1
+        if not all( (i == 2 for i in o.values())):
+            raise Exception("Missing Pair")
+        return ( 'cat', list( (os.path.basename(i), "%s.1.fastq" % i,"%s.2.fastq" % i) for i in o.keys()), 'PE') 
+
+    ### unzipped fastq input
+    fastq_files = walk_dir(base, "*.fastq[12]")
+    if len(fastq_files):
+        o = {}
+        for i in fastq_files:
+            basename = re.sub(r'.fastq[12]$', '', i)
+            o[basename] = o.get(basename, 0) + 1
+        if not all( (i == 2 for i in o.values())):
+            raise Exception("Missing Pair")
+        return ( 'cat', list( (os.path.basename(i), "%s.fastq1" % i,"%s.fastq2" % i) for i in o.keys()), 'PE') 
+
+    ### unzipped txt input
+    fastq_files = walk_dir(base, "*_[12]_sequence.txt")
+    if len(fastq_files):
+        o = {}
+        for i in fastq_files:
+            basename = re.sub(r'_[12]_sequence.txt$', '', i)
+            o[basename] = o.get(basename, 0) + 1
+        if not all( (i == 2 for i in o.values())):
+            raise Exception("Missing Pair")
+        return ( 'cat', list( (os.path.basename(i), "%s_1_sequence.txt" % i,"%s_2_sequence.txt" % i) for i in o.keys()), 'PE') 
+    
+    ### gzipped input
+    fastq_gz_files = walk_dir(base, "*_[12].fastq.gz")
+    if len(fastq_gz_files):
+        o = {}
+        for i in fastq_gz_files:
+            basename = re.sub(r'_[12].fastq.gz$', '', i)
+            o[basename] = o.get(basename, 0) + 1
+        if not all( (i == 2 for i in o.values())):
+            raise Exception("Missing Pair")
+        return ( 'zcat', list( (os.path.basename(i), "%s_1.fastq.gz" % i,"%s_2.fastq.gz" % i) for i in o.keys()), 'PE') 
+
+    ### bzipped input
+    fastq_bz_files = walk_dir(base, "*_[12].fastq.bz")
+    if len(fastq_gz_files):
+        o = {}
+        for i in fastq_gz_files:
+            basename = re.sub(r'_[12].fastq.bz$', '', i)
+            o[basename] = o.get(basename, 0) + 1
+        if not all( (i == 2 for i in o.values())):
+            raise Exception("Missing Pair")
+        return ( 'bzcat', list( (os.path.basename(i), "%s_1.fastq.bz" % i,"%s_2.fastq.bz" % i) for i in o.keys()), 'PE') 
+
+    ### scan for single-end files
+    #############################
+
+    ### unzipped input
+    fastq_files = walk_dir(base, "*.fastq")
+    if len(fastq_files):
+        return ( 'cat', list( (os.path.basename(re.sub(r'.fastq$', '', i)), i) for i in fastq_files), 'SE') 
+
+    ### unzipped input
+    fastq_files = walk_dir(base, "*.fq")
+    if len(fastq_files):
+        return ( 'cat', list( (os.path.basename(re.sub(r'.fq$', '', i)), i) for i in fastq_files), 'SE') 
+
+    ### gzipped input
+    fastq_files = walk_dir(base, "*.fastq.gz")
+    if len(fastq_files):
+        return ( 'zcat', list( (os.path.basename(re.sub(r'.fastq.gz$', '', i)), i) for i in fastq_files), 'SE') 
+
+    ### bzipped input
+    fastq_files = walk_dir(base, "*.fastq.bz")
+    if len(fastq_files):
+        return ( 'bzcat', list( (os.path.basename(re.sub(r'.fastq.bz$', '', i)), i) for i in fastq_files), 'SE') 
+
+    raise Exception("Unable to determine input type")
+
+
+def spreadsheet2dict(spreadFile):
+    """
+    Takes the filename of the spreadsheet, loads the data and organizes
+    it into a dictionary"""
+
+    spreadDict = {}
+    key2field = {}
+    for l, line in enumerate(open(spreadFile)):
+        sl = line.strip().split('\t')
+        if l == 0:
+            for k, key in enumerate(sl):
+                key2field[key] = k
+        else:
+            spreadDict[sl[key2field['analysis_id']]] = sl
+    
+    return (spreadDict, key2field)
+
+
+def spreadsheet2RGdict(spreadFile, analysisID):
+    """Compiles a read group dictionary from the information
+    in the spreadFile for the given analysisID."""
+
+    sD, k2f = spreadsheet2dict(spreadFile)
+
+    try:
+        rec = sD[analysisID]
+    except KeyError:
+        raise Exception('Information for analysis ID %s could not be found in %s' % (analysisID, spreadFile))
+
+    ### build dictionary
+    RG_dict = { 'ID' : '%s:%s' % (rec[k2f['center_name']], analysisID),
+                'CN' : rec[k2f['center_name']],
+                'LB' : 'RNA-Seq:%s:%s' % (rec[k2f['center_name']], rec[k2f['lib_id']]),
+                'PL' : rec[k2f['platform']],
+                'PM' : rec[k2f['platform_model']],
+                'SM' : rec[k2f['specimen_id']],
+                'SI' : rec[k2f['submitted_sample_id']]}
+
+    files = []
+    if 'fastq_files' in k2f:
+        if not rec[k2f['fastq_files']].strip(' ') in ['N/A', 'NA', 'no', '']:
+            files = rec[k2f['fastq_files']].strip(' ').split(' ')
+
+    return (RG_dict, files)
+
+
+def xml2RGdict(xmlfile):
+    
+    ### read xml in
+    root = etree.parse(xmlfile)
+    rtree = root.find('Result')
+
+    ### analysis_id
+    analysis_id = rtree.find('analysis_id').text
+    center = rtree.find('center_name').text
+    try:
+        date_string = rtree.find('analysis_xml/ANALYSIS_SET/ANALYSIS').attrib['analysis_date']
+    except KeyError:
+        date_string = rtree.find('run_xml/RUN_SET/RUN').attrib['run_date']
+    sample_id = rtree.find('sample_id').text
+    submitter_id = rtree.find('legacy_sample_id').text
+    library_id = rtree.find('experiment_xml/EXPERIMENT_SET/EXPERIMENT').attrib['alias']
+    platform = rtree.find('experiment_xml/EXPERIMENT_SET/EXPERIMENT/PLATFORM').getchildren()[0].tag 
+    instrument = rtree.find('experiment_xml/EXPERIMENT_SET/EXPERIMENT/PLATFORM/*/INSTRUMENT_MODEL').text
+
+    ### build dictionary
+    RG_dict = { 'ID' : '%s:%s' % (center, analysis_id),
+                'CN' : center,
+                'DT' : date_string,
+                'LB' : 'RNA-Seq:%s:%s' % (center, library_id),
+                'PL' : platform,
+                'PM' : instrument,
+                'SM' : sample_id,
+                'SI' : submitter_id}
+
+    return RG_dict
+
+
+if __name__ == "__main__":
+
+    parser = argparse.ArgumentParser(description="ICGC RNA-Seq alignment wrapper for STAR alignments.", formatter_class=argparse.ArgumentDefaultsHelpFormatter, usage='%(prog)s [options]', add_help=False)
+    required = parser.add_argument_group("Required input parameters")
+    required.add_argument("--genomeDir", default=None, help="Directory containing the reference genome index", required=True)
+    required.add_argument("--tarFileIn", default=None, help="Input file containing the sequence information", required=True)
+    optional = parser.add_argument_group("optional input parameters")
+    optional.add_argument("--out", default="out.bam", help="Name of the output BAM file")
+    optional.add_argument("--workDir", default="./", help="Work directory")
+    optional.add_argument("--metaDataTab", default=None, help="File containing metadata for the alignment header")
+    optional.add_argument("--analysisID", default=None, help="Analysis ID to be considered in the metadata file")
+    optional.add_argument("--keepJunctions", default=False, action='store_true', help="keeps the junction file as {--out}.junctions")
+    optional.add_argument("--useTMP", default=None, help="environment variable that is used as prefix for temprary data")
+    optional.add_argument("-h", "--help", action='store_true', help="show this help message and exit")
+    star = parser.add_argument_group("STAR input parameters")
+    star.add_argument("--runThreadN", type=int, default=4, help="Number of threads")
+    star.add_argument("--outFilterMultimapScoreRange", type=int, default=1, help="outFilterMultimapScoreRange")
+    star.add_argument("--outFilterMultimapNmax", type=int, default=20, help="outFilterMultimapNmax")
+    star.add_argument("--outFilterMismatchNmax", type=int, default=10, help="outFilterMismatchNmax")
+    star.add_argument("--alignIntronMax", type=int, default=500000, help="alignIntronMax")
+    star.add_argument("--alignMatesGapMax", type=int, default=1000000, help="alignMatesGapMax")
+    star.add_argument("--sjdbScore", type=int, default=2, help="sjdbScore")
+    star.add_argument("--limitBAMsortRAM", type=int, default=0, help="limitBAMsortRAM")
+    star.add_argument("--alignSJDBoverhangMin", type=int, default=1, help="alignSJDBoverhangMin")
+    star.add_argument("--genomeLoad", default="NoSharedMemory", help="genomeLoad")
+    star.add_argument("--genomeFastaFiles", default=None, help="genome sequence in fasta format to rebuild index")
+    star.add_argument("--outFilterMatchNminOverLread", type=float, default=0.33, help="outFilterMatchNminOverLread")
+    star.add_argument("--outFilterScoreMinOverLread", type=float, default=0.33, help="outFilterScoreMinOverLread")
+    star.add_argument("--twopass1readsN", type=int, default=-1, help="twopass1readsN (-1 means all reads are used for remapping)")
+    star.add_argument("--sjdbOverhang", type=int, default=100, help="sjdbOverhang (only necessary for two-pass mode)")
+    star.add_argument("--outSAMstrandField", default="intronMotif", help="outSAMstrandField")
+    star.add_argument("--outSAMattributes", default=["NH", "HI", "NM", "MD", "AS", "XS"], help="outSAMattributes")
+    star.add_argument("--outSAMunmapped", default="Within", help="outSAMunmapped")
+    star.add_argument("--outSAMtype", default=["BAM", "SortedByCoordinate"], help="outSAMtype")
+    star.add_argument("--outSAMheaderHD", default=["@HD", "VN:1.4"], help="outSAMheaderHD")
+    star.add_argument("--outSAMattrRGline", default=None, help="RG attribute line submitted to outSAMattrRGline")
+    star.add_argument("--outSAMattrRGfile", default=None, help="File containing the RG attribute line submitted to outSAMattrRGline")
+    star.add_argument("--outSAMattrRGxml", default=None, help="XML-File in TCGA format to compile RG attribute line")
+    
+    ### check completeness of command line inputs
+    if len(sys.argv) < 2:
+        parser.print_help()
+        sys.exit(0)
+
+    args = parser.parse_args()
+
+    ### some sanity checks on command line parameters
+    if args.metaDataTab is not None:
+        if not os.path.exists(args.metaDataTab):
+            raise Exception("File provided via --metaDataTab does not exist\nFile: %s" % args.metaDataTab)
+        if args.analysisID is None:
+            raise Exception("When providing information in a metadata file, a value for --analysisID is required")
+    if args.outSAMattrRGxml is not None and not os.path.exists(args.outSAMattrRGxml):
+        raise Exception("File provided via --outSAMattrRGxml does not exist\nFile: %s" % args.outSAMattrRGxml)
+    if args.outSAMattrRGfile is not None and not os.path.exists(args.outSAMattrRGfile):
+        raise Exception("File provided via --outSAMattrRGfile does not exist\nFile: %s" % args.outSAMattrRGfile)
+
+    ### handling of input file (unpacking, etc. )
+    if args.useTMP is not None:
+        workdir = tempfile.mkdtemp(dir=os.environ[args.useTMP], prefix="star_inputdir_")
+    else:
+        workdir = tempfile.mkdtemp(dir=args.workDir, prefix="star_inputdir_")
+    if args.tarFileIn.endswith(".gz"):
+        tarcmd = "tar xvzf %s -C %s" % (args.tarFileIn, workdir)
+    elif args.tarFileIn.endswith(".bz"):
+        tarcmd = "tar xvjf %s -C %s" % (args.tarFileIn, workdir)
+    elif args.tarFileIn.endswith(".tar"):
+        tarcmd = "tar xvf %s -C %s" % (args.tarFileIn, workdir)
+    elif args.tarFileIn.endswith(".sra"):
+        tarcmd = "fastq-dump --gzip --split-3 --outdir %s %s" % (workdir, args.tarFileIn)
+    else:
+        raise Exception("Unknown input file extension for file %s" % (args.tarFileIn))
+    subprocess.check_call(tarcmd, shell=True)
+    
+    ### collect fastq information from extraction dir
+    align_sets = scan_workdir(os.path.abspath(workdir))
+    
+    ### process read group information
+    files = []
+    if args.metaDataTab is not None:
+        (RG_dict, files_tmp) = spreadsheet2RGdict(args.metaDataTab, args.analysisID) 
+        files.extend(files_tmp)
+    elif args.outSAMattrRGxml is not None:
+        RG_dict = xml2RGdict(args.outSAMattrRGxml)
+    elif args.outSAMattrRGline is not None:
+        RG_dict = dict([(x.split(':', 1)[0], x.split(':', 1)[1]) for x in args.outSAMattrRGline.split()])
+    elif args.outSAMattrRGfile is not None:
+        _fh = open(args.outSAMattrRGfile, 'r')
+        RG_dict = dict([(x.split(':', 1)[0], x.split(':', 1)[1]) for x in _fh.next().strip().split()])
+        _fh.close()
+    else:
+        RG_dict = {'ID' : '', 'SM' : ''}
+
+    ### post-process RG-dict to comply with STAR conventions
+    for key in RG_dict:
+        sl = RG_dict[key].split(' ')
+        if len(sl) > 1:
+            RG_dict[key] = '"%s"' % RG_dict[key]
+
+    ### use list of fastq input files for whitelisting
+    if len(files) > 0:
+        align_sets = (align_sets[0], [x for x in align_sets[1] if (re.sub('(_[12]){,1}.fastq(.(gz|bz2|bz))*', '', os.path.basename(x[1])) in files)], align_sets[2])
+        if len(align_sets[1]) == 0:
+            print >> sys.stderr, 'All input files have been filtered out - no input remaining. Terminating.'
+            sys.exit()
+
+    ### use filename stub as read group label
+    RG_dict['RG'] = []
+    for fn in [x[1] for x in align_sets[1]]:
+        fn_stub = re.sub('(_[12]){,1}.fastq(.(gz|bz2|bz))*', '', os.path.basename(fn))
+        fn_stub = re.sub('(_[12]){,1}_sequence.txt(.(gz|bz2|bz))*', '', fn_stub)
+        fn_stub = re.sub('_read[12]', '', fn_stub)
+        fn_stub = re.sub('_R[12]_001$', '', fn_stub)
+        RG_dict['RG'].append(fn_stub)
+
+    ### prepare template string
+    if align_sets[2] == 'PE':
+        read_str = '${fastq_left} ${fastq_right}'
+    else:
+        read_str = '${fastq_left}'
+    
+    ### simulate two pass alignment until STAR fully implements this
+    if args.twopass1readsN != 0:
+        
+        ### run first round of alignments and only record junctions
+        align_template_str_1st = """STAR \
+--genomeDir ${genomeDir} --readFilesIn %s \
+--runThreadN ${runThreadN} \
+--outFilterMultimapScoreRange ${outFilterMultimapScoreRange} \
+--outFilterMultimapNmax ${outFilterMultimapNmax} \
+--outFilterMismatchNmax ${outFilterMismatchNmax} \
+--alignIntronMax ${alignIntronMax} \
+--alignMatesGapMax ${alignMatesGapMax} \
+--sjdbScore ${sjdbScore} \
+--alignSJDBoverhangMin ${alignSJDBoverhangMin} \
+--genomeLoad ${genomeLoad} \
+--readFilesCommand ${readFilesCommand} \
+--outFilterMatchNminOverLread ${outFilterMatchNminOverLread} \
+--outFilterScoreMinOverLread ${outFilterScoreMinOverLread} \
+--sjdbOverhang ${sjdbOverhang} \
+--outSAMstrandField ${outSAMstrandField} \
+--outSAMtype None \
+--outSAMmode None""" % read_str
+
+        if args.twopass1readsN > 0:
+            align_template_str_1st += " --readMapNumber %i" % args.twopass1readsN
+
+        cmd = string.Template(align_template_str_1st).safe_substitute({
+            'genomeDir' : os.path.abspath(args.genomeDir),
+            'runThreadN' : args.runThreadN,
+            'outFilterMultimapScoreRange' : args.outFilterMultimapScoreRange,
+            'outFilterMultimapNmax' : args.outFilterMultimapNmax,
+            'outFilterMismatchNmax' : args.outFilterMismatchNmax,
+            'fastq_left' : ','.join([os.path.join(x[0], x[1]) for x in align_sets[1]]),
+            'alignIntronMax' : args.alignIntronMax,
+            'alignMatesGapMax': args.alignMatesGapMax,
+            'sjdbScore': args.sjdbScore,
+            'alignSJDBoverhangMin' : args.alignSJDBoverhangMin,
+            'genomeLoad' : args.genomeLoad,
+            'readFilesCommand' : align_sets[0],
+            'outFilterMatchNminOverLread' : args.outFilterMatchNminOverLread,
+            'outFilterScoreMinOverLread' : args.outFilterScoreMinOverLread,
+            'sjdbOverhang' : args.sjdbOverhang,
+            'outSAMstrandField' : args.outSAMstrandField
+        })
+        if align_sets[2] == 'PE':
+            cmd = string.Template(cmd).substitute({
+                'fastq_right' : ','.join([os.path.join(x[0], x[2]) for x in align_sets[1]])
+            })
+
+        ### take temp directory from environment variable
+        if args.useTMP is not None:
+            align_dir_1st = os.path.abspath( tempfile.mkdtemp(dir=os.environ[args.useTMP], prefix="star_aligndir_1st_") )
+            genome_dir_1st = os.path.abspath( tempfile.mkdtemp(dir=os.environ[args.useTMP], prefix="star_genomedir_1st_") )
+        else:
+            align_dir_1st = os.path.abspath( tempfile.mkdtemp(dir=args.workDir, prefix="star_aligndir_1st_") )
+            genome_dir_1st = os.path.abspath( tempfile.mkdtemp(dir=args.workDir, prefix="star_genomedir_1st_") )
+        print "Running", cmd
+        subprocess.check_call(cmd, shell=True, cwd=align_dir_1st)
+
+        ### build index using provided genome fasta as well as junctions from first run
+        cmd = """STAR --runMode genomeGenerate --genomeDir %s \
+--genomeFastaFiles %s \
+--sjdbOverhang %i \
+--runThreadN %i \
+--sjdbFileChrStartEnd %s""" % (genome_dir_1st, args.genomeFastaFiles, args.sjdbOverhang, args.runThreadN, os.path.join(align_dir_1st, 'SJ.out.tab')) 
+        print "Running", cmd
+        subprocess.check_call(cmd, shell=True, cwd=align_dir_1st)
+
+        ### replace index for the second run with the one currently built
+        genome_dir = genome_dir_1st
+    else:
+        genome_dir = os.path.abspath(args.genomeDir)
+
+
+    align_template_str = """STAR \
+--genomeDir ${genomeDir} --readFilesIn %s \
+--runThreadN ${runThreadN} \
+--outFilterMultimapScoreRange ${outFilterMultimapScoreRange} \
+--outFilterMultimapNmax ${outFilterMultimapNmax} \
+--outFilterMismatchNmax ${outFilterMismatchNmax} \
+--alignIntronMax ${alignIntronMax} \
+--alignMatesGapMax ${alignMatesGapMax} \
+--sjdbScore ${sjdbScore} \
+--alignSJDBoverhangMin ${alignSJDBoverhangMin} \
+--genomeLoad ${genomeLoad} \
+--limitBAMsortRAM ${limitBAMsortRAM} \
+--readFilesCommand ${readFilesCommand} \
+--outFilterMatchNminOverLread ${outFilterMatchNminOverLread} \
+--outFilterScoreMinOverLread ${outFilterScoreMinOverLread} \
+--sjdbOverhang ${sjdbOverhang} \
+--outSAMstrandField ${outSAMstrandField} \
+--outSAMattributes ${outSAMattributes} \
+--outSAMunmapped ${outSAMunmapped} \
+--outSAMtype ${outSAMtype} \
+--outSAMheaderHD ${outSAMheaderHD}""" % read_str
+
+#--twopass1readsN ${twopass1readsN} \
+
+    cmd = string.Template(align_template_str).safe_substitute({
+        'genomeDir' : genome_dir,
+        'runThreadN' : args.runThreadN,
+        'fastq_left' : ','.join([os.path.join(x[0], x[1]) for x in align_sets[1]]), #os.path.abspath(pair[1]),
+        'outFilterMultimapScoreRange' : args.outFilterMultimapScoreRange,
+        'outFilterMultimapNmax' : args.outFilterMultimapNmax,
+        'outFilterMismatchNmax' : args.outFilterMismatchNmax,
+        'alignIntronMax' : args.alignIntronMax,
+        'alignMatesGapMax': args.alignMatesGapMax,
+        'sjdbScore': args.sjdbScore,
+        'alignSJDBoverhangMin' : args.alignSJDBoverhangMin,
+        'genomeLoad' : args.genomeLoad,
+        'limitBAMsortRAM' : args.limitBAMsortRAM,
+        'readFilesCommand' : align_sets[0],
+        'outFilterMatchNminOverLread' : args.outFilterMatchNminOverLread,
+        'outFilterScoreMinOverLread' : args.outFilterScoreMinOverLread,
+        'sjdbOverhang' : args.sjdbOverhang,
+        'outSAMstrandField' : args.outSAMstrandField,
+        'outSAMattributes' : " ".join(args.outSAMattributes),
+        'outSAMunmapped' : args.outSAMunmapped, 
+        'outSAMtype' : " ".join(args.outSAMtype),
+        'outSAMheaderHD' : " ".join(args.outSAMheaderHD)
+    })
+#        'twopass1readsN' : args.twopass1readsN,
+    if align_sets[2] == 'PE':
+        cmd = string.Template(cmd).substitute({
+            'fastq_right' : ','.join([os.path.join(x[0], x[2]) for x in align_sets[1]]) # os.path.abspath(pair[2]),
+        })
+
+    ### convert RG_dict into formatted RG line
+    RG_line = []
+    for r, readgroup in enumerate(align_sets[1]):
+        if 'RG' in RG_dict:
+            tmp = 'ID:%s:%s' % (RG_dict['ID'], RG_dict['RG'][r])
+        else:
+            tmp = 'ID:%s:%s' % (RG_dict['ID'], readgroup[0])
+        if len(RG_dict) > 1:
+            tmp += '\t'
+            tmp += '\t'.join(['%s:%s' % (key, RG_dict[key]) for key in RG_dict if key not in ['ID', 'RG', 'SI']])
+        ### add read group label
+        if 'RG' in RG_dict and 'CN' in RG_dict:
+            tmp += '\tPU:%s:%s' % (RG_dict['CN'], RG_dict['RG'][r])
+        RG_line.append('%s' % tmp)
+    cmd += ' --outSAMattrRGline %s' % ' , '.join(RG_line)
+
+    ### handle comment lines
+    comment_file = None
+    if 'SI' in RG_dict:
+        if args.useTMP is not None:
+            comment_file = os.path.abspath( tempfile.mkstemp(dir=os.environ[args.useTMP], prefix="star_comments_")[1] )
+        else:
+            comment_file = os.path.abspath( tempfile.mkstemp(dir=args.workDir, prefix="star_comments_")[1] )
+        
+        fd_com = open(comment_file, 'w')
+        fd_com.write('@CO\tsubmitter_sample_id:%s\n' % RG_dict['SI'])
+
+        fd_com.flush()
+        fd_com.close()
+        
+        cmd += ' --outSAMheaderCommentFile %s' % comment_file
+
+
+    ### take temp directory from environment variable
+    if args.useTMP is not None:
+        align_dir = os.path.abspath( tempfile.mkdtemp(dir=os.environ[args.useTMP], prefix="star_aligndir_") )
+    else:
+        align_dir = os.path.abspath( tempfile.mkdtemp(dir=args.workDir, prefix="star_aligndir_") )
+    print "Running", cmd
+    subprocess.check_call(cmd, shell=True, cwd=align_dir)
+
+    ### move output file
+    if 'BAM' in args.outSAMtype and 'SortedByCoordinate' in args.outSAMtype:
+        shutil.move(os.path.join(align_dir, 'Aligned.sortedByCoord.out.bam'), args.out)
+    elif 'BAM' in args.outSAMtype and 'Unsorted' in args.outSAMtype:
+        shutil.move(os.path.join(align_dir, 'Aligned.out.bam'), args.out)
+    else:
+        raise Exception('STAR output file could not be determined') 
+
+    ### move junctions if to be kept
+    if args.keepJunctions:
+        shutil.move(os.path.join(align_dir, 'SJ.out.tab'), args.out + '.junctions')
+
+    ### clean up working directory
+    shutil.rmtree(workdir)
+    shutil.rmtree(align_dir)
+    if args.twopass1readsN != 0:
+        shutil.rmtree(align_dir_1st)
+        shutil.rmtree(genome_dir_1st)