FindDNAFusion/generateBAM.pl at main · jml-bioinfo/FindDNAFusion · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
#!/usr/local/bin/perl


use strict;
use Getopt::Long;
use Cwd qw(cwd);
use POSIX qw(strftime);


=head1 NAME

runGeneFuse.pl


=head2 SYNOPSIS

perl generateBAM.pl -i input-dir -o output-dir


=head3 DESCRIPTION

This PERL script to generate BAM files for multiple samples.


=head4 AUTHOR

Xiaokang Pan (Xiaokang.Pan@osumc.edu)


=head5 LAST UPDATE

01/29/2023

=cut


my $usage = <<EOS;
   Usage: perl $0 -i inpath -r ref_genome [-options]

   -i|inpath	  [string] (input sequence directory storing fastq files)
   -r|ref_genome  [string] (reference genome in FASTA format)
   -o|outpath	  [string] (output directory)
   -c|cpus        [string] CPU number per sample to be used
   -h|help        (help information)

EOS

my ($inpath, $outpath, $ref_genome, $cpus, $help);
GetOptions (
  "inpath=s"     => \$inpath,        # input dir
  "ref_genome=s" => \$ref_genome,    # reference genome
  "outpath:s"    => \$outpath,       # output dir
  "cpus:s"       => \$cpus,          # number of CPUs
  "help:s"       => \$help	     # help information
);
die "\n$usage\n" if ($help or
                     !defined($inpath) or
                     !defined($ref_genome));

my $current_path = cwd;

if ($inpath !~ /^\/.+\/.+$/) {
    $inpath =  $current_path."/".$inpath;
}

if (!$outpath) {
    system("mkdir mbam");
    $outpath = $current_path."/mbam";

} elsif ($outpath !~ /^\/.+\/.+$/) {
    $outpath = $current_path."/".$outpath;
}

if (!$cpus) {
    $cpus = 16;
}

my @files = <$inpath/*.fastq.gz>;

my %samples;
foreach my $file (@files) {

    my $sample_name;

    if ($file =~ /^.+\/(.+)_R\d\.fastq\.gz$/) {
        $sample_name = $1;
        $samples{$sample_name} = 1;
    }
}

#my $ref_file = "database/hg19.fasta";

# run BWA MEM to map on Human Genome
foreach my $sample (sort keys %samples) {

   next if ($sample =~ /-RNA$/ || $sample =~ /-R$/);

   my $r1_fastq_file = "$inpath/$sample"."_R1.fastq.gz";
   my $r2_fastq_file = "$inpath/$sample"."_R2.fastq.gz";
   system("bwa mem -M -t $cpus $ref_genome $r1_fastq_file $r2_fastq_file | samtools sort -@ $cpus -o $outpath/$sample.bam - &");
}
system("sleep 60m");

# check if mapping completed.
# if not, wait for a few minutes
my @tmp_files = <$outpath/*.tmp.*>;
while (@tmp_files) {
    system("sleep 5m");
    @tmp_files = <$outpath/*.tmp.*>;
}

# make index
my @bam_files = <$outpath/*.bam>;
foreach my $bam_file (@bam_files) {

    system("samtools index $bam_file &");
}
system("sleep 10m");
print "BAM generation is completed.\n\n";