Background information
I've been asked to write a little Perl script that allows genomic data to be screened against reference files in order to determine locations of specific mutations.
The input file format look like this (information are tab-separated). In the script it's referred to as Funestus_NON_SYN.txt.
KB669306 9962 . C T 520.77 PASS AC=13;AF=0.929;AN=14;BaseQRankSum=0.540;DP=103;Dels=0.00;EFF=SYNONYMOUS_CODING(LOW|SILENT|cgC/cgT|R15||AFUN012013|||AFUN012013-RA|1|1);FS=0.000;HaplotypeScore=0.0000;MQ0=0;MQRankSum=0.231;ReadPosRankSum=0.694;set=variant-variant2-variant3-variant5-variant6-variant7-variant10 GT:AD:DP:GQ:PL 1/1:0,18:18:42:549,42,0 1/1:0,12:12:33:399,33,0 1/1:0,13:13:39:505,39,0 ./. 0/1:9,4:13:99:104,0,297 1/1:0,17:17:45:570,45,0 1/1:0,18:18:51:624,51,0 ./. ./. 1/1:0,12:12:33:420,33,0
KB669306 10439 . C T 668.77 PASS AC=15;AF=0.938;AN=16;BaseQRankSum=-0.643;DP=139;Dels=0.00;EFF=SYNONYMOUS_CODING(LOW|SILENT|gcC/gcT|A174||AFUN012013|||AFUN012013-RA|1|1);FS=0.000;MQ0=0;MQRankSum=-0.746;ReadPosRankSum=-1.106;set=variant-variant3-variant5-variant6-variant7-variant8-variant9-variant10 GT:AD:DP:GQ:PL 1/1:0,22:22:51:697,51,0 ./. 1/1:0,20:20:48:667,48,0 ./. 0/1:14,12:25:99:317,0,403 1/1:0,15:15:36:469,36,0 1/1:0,16:16:36:499,36,0 1/1:0,13:13:30:415,30,0 1/1:0,14:14:30:412,30,0 1/1:0,13:13:36:492,36,0
KB668289 903 . T C 577.77 PASS AC=14;AF=0.875;AN=16;DP=173;Dels=0.00;EFF=DOWNSTREAM(MODIFIER||3412|||AFUN005870|||AFUN005870-RA||1),SYNONYMOUS_CODING(LOW|SILENT|ctA/ctG|L578||AFUN005868|||AFUN005868-RA|4|1),UPSTREAM(MODIFIER||2667|||AFUN005869|||AFUN005869-RA||1);MQ0=0;set=variant-variant2-variant3-variant5-variant6-variant7-variant9-variant10 GT:AD:DP:GQ:PL 1/1:0,20:20:45:606,45,0 1/1:0,19:19:42:588,42,0 1/1:0,22:22:54:720,54,0 ./. 1/1:0,29:29:63:882,63,0 0/1:16,13:28:99:401,0,391 1/1:0,13:13:30:425,30,0./. 1/1:0,26:26:57:785,57,0 0/1:8,7:15:99:223,0,243
KB668289 1224 . A C 572.77 PASS AC=10;AF=0.833;AN=12;DP=122;Dels=0.00;EFF=DOWNSTREAM(MODIFIER||3091|||AFUN005870|||AFUN005870-RA||1),SYNONYMOUS_CODING(LOW|SILENT|cgT/cgG|R471||AFUN005868|||AFUN005868-RA|4|1),UPSTREAM(MODIFIER||2346|||AFUN005869|||AFUN005869-RA||1);HaplotypeScore=0.0000;MQ0=0;set=variant-variant3-variant6-variant7-variant8-variant9 GT:AD:DP:GQ:PL 1/1:1,20:20:42:601,42,0 ./. 1/1:0,26:26:63:865,63,0 ./. ./. 0/1:9,10:19:99:276,0,236 1/1:0,17:17:42:583,42,0 0/1:14,5:19:99:119,0,406 1/1:0,20:20:45:630,45,0 ./.
Now, I have to screen the first column against files containing genes in the following format. Referred to as gsts-funestus.gff in the script
# start gene FUNEGST002
KB668289 AUGUSTUS gene 410926 411627 0.18 + . FUNEGST002; gene_name "GSTd3"; auto "g17128_modified";
KB668289 AUGUSTUS transcript 410926 411627 0.18 + . FUNEGST002.t1
KB668289 AUGUSTUS start_codon 410926 410928 . + 0 transcript_id "FUNEGST002.t1"; gene_id "FUNEGST002";
KB668289 AUGUSTUS intron 411058 411126 0.99 + . transcript_id "FUNEGST002.t1"; gene_id "FUNEGST002";
KB668289 AUGUSTUS CDS 410926 411057 0.58 + 0 transcript_id "FUNEGST002.t1"; gene_id "FUNEGST002";
KB668289 AUGUSTUS CDS 411127 411627 0.74 + 0 transcript_id "FUNEGST002.t1"; gene_id "FUNEGST002";
KB668289 AUGUSTUS stop_codon 411625 411627 . + 0 transcript_id "FUNEGST002.t1"; gene_id "FUNEGST002";
# end gene FUNEGST002
Basically I have to extract the first line containing gene
in the 3rd column and print a hybrid output in the following format.
GSTd3 KB668289 903 T=>C AC=14;AF=0.875;AN=16;DP=173;Dels=0.00;EFF=DOWNSTREAM(MODIFIER||3412|||AFUN005870|||AFUN005870-RA||1),SYNONYMOUS_CODING(LOW|SILENT|ctA/ctG|L578||AFUN005868|||AFUN005868-RA|4|1),UPSTREAM(MODIFIER||2667|||AFUN005869|||AFUN005869-RA||1);MQ0=0;set=variant-variant2-variant3-variant5-variant6-variant7-variant9-variant10 GT:AD:DP:GQ:PL 1/1:0,20:20:45:606,45,0 1/1:0,19:19:42:588,42,0 1/1:0,22:22:54:720,54,0 ./. 1/1:0,29:29:63:882,63,0 0/1:16,13:28:99:401,0,391 1/1:0,13:13:30:425,30,0./. 1/1:0,26:26:57:785,57,0 0/1:8,7:15:99:223,0,243
GSTd3 KB668289 1224 A=>C AC=10;AF=0.833;AN=12;DP=122;Dels=0.00;EFF=DOWNSTREAM(MODIFIER||3091|||AFUN005870|||AFUN005870-RA||1),SYNONYMOUS_CODING(LOW|SILENT|cgT/cgG|R471||AFUN005868|||AFUN005868-RA|4|1),UPSTREAM(MODIFIER||2346|||AFUN005869|||AFUN005869-RA||1);HaplotypeScore=0.0000;MQ0=0;set=variant-variant3-variant6-variant7-variant8-variant9 GT:AD:DP:GQ:PL 1/1:1,20:20:42:601,42,0 ./. 1/1:0,26:26:63:865,63,0 ./. ./. 0/1:9,10:19:99:276,0,236 1/1:0,17:17:42:583,42,0 0/1:14,5:19:99:119,0,406 1/1:0,20:20:45:630,45,0 ./.
The Script
My script contains this part to read and process all lines in the input file.
#!/usr/bin/perl -w
use strict;
use warnings;
print "Enter your input file containing all SNPs [e.g. Funestus_NON_SYN.txt]:\t";
chomp (my $input_file = <STDIN>);
# Open and import the information found in the files #
my @input = &InputFileHandler($input_file);
# Store the gene information in a hash with the key as scaffold and location #
my %genes = ();
foreach (@input)
{
next if $_ !~ /^KB.*/;
my($scaffold,$location,$dot,$start,$end,$unknown,$pass,$rest) = split(/\t/, $_, 8);
my $mutation = "$start=>$end";
my $key = $scaffold."_".$location;
my $value = "$scaffold\t$location\t$mutation\t$rest";
$genes{$key} = $value;
}
The next part allows the user to define the number of reference files and screens out input file against it. Furthermore, it prints the results in a tab-separated output file.
# Open the output file
my $output_file = &GenerateFilename();
open my $fh_out,">",$output_file.".tsv" or die "Cannot open $output_file: $!\n";
# Print the header in the CSV file
print $fh_out "Gene\tSNP\tLocation\tMutation\tRemaining information\n###\n";
print "Enter the number of reference files [e.g. gff files]:\t";
chomp (my $i = <STDIN>);
# Generic loop that let's the user enter as many reference files as required #
for (my $j = 1; $j <= $i;$j++)
{
print "Enter filename $j [e.g. gsts-funestus.gff]:\t";
chomp (my $file = <STDIN>);
my @file_raw = &InputFileHandler($file);
my @file_genes = &ScanArray(@file_raw);
my @file_genes_print = &ScanPrint(@file_genes);
print $fh_out "$_\n" foreach (@file_genes_print);
print $fh_out "\n###\n";
}
Finally, my subroutines are the following.
sub GenerateFilename # Generates a filename based on the current date and time.
{
my @timeData = localtime(time);
my $filename = join('_', @timeData);
return $filename;
}
sub InputFileHandler # Opens a file and saves its content in an array
{
my $file = shift;
my @lines;
open my $fh,"<",$file or die "Cannot open $file: $!\n";
while (my $line = <$fh>)
{
chomp($line);
push (@lines,$line);
}
close ($fh);
return @lines;
}
sub Information # Separates the mutation from the remaining information
{
my $rest = shift;
my @rest = split (/\t/,$rest);
my $aa_change = "$rest[1]=>$rest[2]";
my $return = "$aa_change,$rest[5]";
return $return;
}
sub ScanArray # Scans a file array for lines containing SNPs, discards rest
{
my @array = @_;
my @return;
foreach (@array)
{
next if $_ !~ /^KB.*\t.*\t[Gg][Ee][Nn][Ee]/;
my ($scaffold,$second,$type,$start,$end,$sixth,$seventh,$eigth,$info) = split (/\t/,$_);
my ($unknown,$gene_name,$auto) = split (/; /,$info,3);
$gene_name =~ s/.*"(.*?)".*/$1/s;
my $gene = $scaffold."_".$gene_name;
push (@return,$gene);
}
return @return;
}
sub ScanPrint # Filters for genes that contain SNPs
{
my @array = @_;
my @return;
foreach (@array)
{
my ($scaffold,$gene_name) = split (/_/,$_,2);
foreach my $key (keys %genes)
{
next if $key !~ /$scaffold/;
my $line = "$gene_name\t$genes{$key}";
push (@return, $line);
}
}
return @return;
}
Question
Do you people have any suggestions or ways of improving this script. It all works and no warnings have been returned with all search queries I've done (I used different files of the same format).
EDIT 1
As pointed out to me by @Edward: Yes, my output should be the two lines. I have adjusted it in the main part of the question.