<commit_before>#!/usr/bin/env python
# -*- coding: utf-8 -*-

from subprocess import call
import os

import argparse
parser = argparse.ArgumentParser()

parser.add_argument("-i", "--input", help="BAM file (can be the location on S3)")


args = parser.parse_args()
bam_file = args.input


print(bam_file)

#mock up
# bam_file = "../../../input/bam/test.recal.bam"


input_folder = "/home/ubuntu/projects/input/bam"

base=os.path.basename(bam_file)

if bam_file.startswith('s3://'):
    #download file to input folder
    command = "s3cmd get %s %s/" % (bam_file, input_folder)
    output = call(command, shell=True)
    print(output)
    bam_file = "%s/%s" % (input_folder, base)

# base_name = os.path.splitext(base)[0]
base_name = base.split('.')[0]
print(base_name)

memory_use = "15g"

human_reference = "/home/ubuntu/projects/input/b37/human_g1k_v37.fasta" #84 features
human_reference = "/home/ubuntu/projects/input/grch37/d5/hs37d5.fa" #86 features

#create one folder per sample
output_folder = "/home/ubuntu/projects/output/reports/bam/%s" % (base_name)
if not os.path.exists(output_folder):
    os.makedirs(output_folder)

samtools_dir = "/home/ubuntu/projects/programs/samtools-1.3.1"
picard_dir = "/home/ubuntu/projects/programs/picard"
gatk_dir = "/home/ubuntu/projects/programs/gatk"
qualimap_dir = "/home/ubuntu/projects/programs/qualimap/qualimap_v2.2"

#s3
#if bam file start with s3 download from s3

if not os.path.exists(bam_file+'.bai'):
    print('Indexing BAM')
    #index bam
    command = "%s/samtools index %s" % (samtools_dir, bam_file)
    output = call(command, shell=True)
    print(output)

#fastqc
#done already!

print('Running featureCounts')
#featureCounts
command = """/home/ubuntu/projects/programs/subread-1.5.1-Linux-x86_64/bin/featureCounts --donotsort -T 4 -p \
-a /home/ubuntu/projects/input/gtf/Homo_sapiens.GRCh37.75.gtf \
-o %s/%s.featureCounts.txt \
%s""" % (output_folder, base_name, bam_file)
# output = call(command, shell=True)
# print(output)

#bamtools
print('Running bamtools')
command = """/home/ubuntu/projects/programs/bamtools/bin/bamtools stats -in %s > %s/%s.bamtools.stats.txt""" % (bam_file, output_folder, base_name)
# output = call(command, shell=True)
# print(output)

print('Running DepthOfCoverage')
#gatk DepthOfCoverage
command = """
java -Xmx15g -jar %s/GenomeAnalysisTK.jar -T DepthOfCoverage \
-I %s \
-R %s \
-o %s/%s.DepthOfCoverage.txt \
-nt 4 \
--omitIntervalStatistics \
-ct 5 -ct 10 -ct 20 -ct 30 \
-log %s/%s.DepthofCoverage.log \
""" % (gatk_dir, bam_file, human_reference, output_folder, base_name, output_folder, base_name)
output = call(command, shell=True)
print(output)

#qualimap BamQC
print('Running qualimap BamQC')
command = """%s/qualimap bamqc \
--java-mem-size=15G \
-bam %s \
-outdir %s \
-nt 4
""" % (qualimap_dir, bam_file, output_folder)
output = call(command, shell=True)
print(output)
 
#samtools flagstat
print('Running samtools flagstat')
command = """%s/samtools flagstat %s > %s/%s.samtools.flagstat.txt
""" % (samtools_dir, bam_file, output_folder, base_name)<commit_msg>bam genome stats<commit_after>#!/usr/bin/env python
# -*- coding: utf-8 -*-

from subprocess import call
import os

import argparse
parser = argparse.ArgumentParser()

parser.add_argument("-i", "--input", help="BAM file (can be the location on S3)")


args = parser.parse_args()
bam_file = args.input


print(bam_file)

#mock up
# bam_file = "../../../input/bam/test.recal.bam"


input_folder = "/home/ubuntu/projects/input/bam"

base=os.path.basename(bam_file)

if bam_file.startswith('s3://'):
    #download file to input folder
    command = "s3cmd get %s %s/" % (bam_file, input_folder)
    output = call(command, shell=True)
    print(output)
    bam_file = "%s/%s" % (input_folder, base)

# base_name = os.path.splitext(base)[0]
base_name = base.split('.')[0]
print(base_name)

memory_use = "15g"

human_reference = "/home/ubuntu/projects/input/b37/human_g1k_v37.fasta" #84 features
human_reference = "/home/ubuntu/projects/input/grch37/d5/hs37d5.fa" #86 features

#create one folder per sample
output_folder = "/home/ubuntu/projects/output/reports/bam/%s" % (base_name)
if not os.path.exists(output_folder):
    os.makedirs(output_folder)

samtools_dir = "/home/ubuntu/projects/programs/samtools-1.3.1"
picard_dir = "/home/ubuntu/projects/programs/picard"
gatk_dir = "/home/ubuntu/projects/programs/gatk"
qualimap_dir = "/home/ubuntu/projects/programs/qualimap/qualimap_v2.2"

#s3
#if bam file start with s3 download from s3

if not os.path.exists(bam_file+'.bai'):
    print('Indexing BAM')
    #index bam
    command = "%s/samtools index %s" % (samtools_dir, bam_file)
    output = call(command, shell=True)
    print(output)

#fastqc
#done already!

print('Running featureCounts')
#featureCounts
command = """/home/ubuntu/projects/programs/subread-1.5.1-Linux-x86_64/bin/featureCounts --donotsort -T 4 -p \
-a /home/ubuntu/projects/input/gtf/Homo_sapiens.GRCh37.75.gtf \
-o %s/%s.featureCounts.txt \
%s""" % (output_folder, base_name, bam_file)
# output = call(command, shell=True)
# print(output)

#bamtools
print('Running bamtools')
command = """/home/ubuntu/projects/programs/bamtools/bin/bamtools stats -in %s > %s/%s.bamtools.stats.txt""" % (bam_file, output_folder, base_name)
# output = call(command, shell=True)
# print(output)

print('Running DepthOfCoverage')
#gatk DepthOfCoverage
command = """
java -Xmx15g -jar %s/GenomeAnalysisTK.jar -T DepthOfCoverage \
-I %s \
-R %s \
-o %s/%s.DepthOfCoverage.txt \
-nt 4 \
--omitIntervalStatistics \
-ct 5 -ct 10 -ct 20 -ct 30 \
-log %s/%s.DepthofCoverage.log \
""" % (gatk_dir, bam_file, human_reference, output_folder, base_name, output_folder, base_name)
# output = call(command, shell=True)
# print(output)

#qualimap BamQC
print('Running qualimap BamQC')
command = """%s/qualimap bamqc \
--java-mem-size=15G \
-bam %s \
-outdir %s \
-nt 4
""" % (qualimap_dir, bam_file, output_folder)
output = call(command, shell=True)
print(output)
 
#samtools flagstat
print('Running samtools flagstat')
command = """%s/samtools flagstat %s > %s/%s.samtools.flagstat.txt
""" % (samtools_dir, bam_file, output_folder, base_name)
{'fraction_non_alphanumeric': 0.16029913835148757, 'fraction_numerical': 0.01820842139489514, 'mean_word_length': 4.429832303618712, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1031740', 'n_tokens_mistral': 2511, 'n_tokens_neox': 2386, 'n_words': 475}