Help In Writing A Perl Script To Call Somatic Variants

Question

[shilpy]

The script should take 2 input vcf files(normal & tumor) and a output file.If a specific variant is found at a position in the tumor sample vcf file but there is no corresponding variant in the normal sample vcf the results will be written to the output file. Please help or guide for a clue how to work on this.

[ff.cc.cc]

Hi. I suggest you to involve a bioinformatician in your project and then ask biostar for tips, not for the whole work...

[Sean Davis]

Just a note that what you are describing is not likely to be a very effective somatic variant caller. There are MANY tools for somatic variant calling and none of them use the approach you are suggesting, for pretty good reason.

[shilpy]

It may not be very effective way but still it is used as I am asked to do in my class.

[Sean Davis]

There are often many approaches to answer a question. "As long as we don't ever evaluate the answers critically, all methods are equivalent as long as they generate results." Many of us experience this sentiment--resist it!

[Michael]

I think it is also problematic that he wants us to do his homework.

[Obi Griffith]

Agreed. It is generally poor practice to just post your homework on a forum. But, at least the poster has admitted it and made some attempt to solve. Although, as it is, the question is more about very beginner Perl scripting and not really a bioinformatic question. I wonder if you wouldn't have more luck reading some forums and tutorials specifically on Perl. Once you have figured out how to make a perl script, understand conditional statements, hash data structures, etc come back for help with the components of this problem that are really bioinformatics (e.g., understanding VCF format).

[shilpy]

Hello guys thanks for your comments, finally I got my perl script working, and it is working fine. Now I am able to compare 2 files and get somatic variants. I came to this forum thinking people will help me, unfortunately I got few tips only and much of discouragement. I want to thank members who commented sarcastically too, as your negative comments helped me to work more hard and get desired output.

[Michael]

I think the correct guidance for getting a clue how to work is to google for "vcf diff", or "vcf variant" or similar, or beginning by studying the available documentation of the vcf format and the available tools, e.g. vcf-tools. Even thought about reading wikipedia articles? In conclusion, the right way to go is to improve ones general web/tech literacy first (learn to search), to make optimal use of the available resources.

[Jordan]

I agree with Michael. Must of the things you need, are usually available online. If the work you are doing is part of a HW, you will need to spend a bit more time to understand things.

Answer 1

[Michael]

This can be done by vcf-compare and vcf-isec.

Answer 2

[shilpy]

!/usr/bin/perl
use strict;
use warnings;

my $sift_string;
my $sift_string2;

open VCF, "./normal_variant_filtered.vcf" or die "Unable to open file: $!";
open VCF1,"./tumor_variant_filtered.vcf" or die "Unable to open file:$!";
my @lines = <VCF>;
print "@lines\n";
my @lines2 = <VCF1>;
print "@lines\n";

if ("normal_variant_filtered.vcf" eq "tumor_variant_filtered.vcf")  {
    print "They're equal\n";
}else{
    print "They are not equal\n";
}

[Obi Griffith]

I think your script needs some work. You will need to go through each line of the tumor vcf file (and perhaps create a hash of the chromosome, start and end position of each variant) and compare to each line in the normal vcf file. Currently your conditional statement is just checking to see if one string is equivalent to another different string (it can't be).

[shilpy]

Please help me how to go about hash part and comparision.

[shilpy]

my $sift_string; my $sift_string2;

open VCF, "./normal_variant_filtered.vcf" or die "Unable to open file: $!"; open VCF1,"./tumor_variant_filtered.vcf" or die "Unable to open file:$!"; my @lines = <vcf>; my @lines1= splice @lines,0,27; #print "@lines\n"; print"**********************************************\n"; my @lines2 = <vcf1>; my @lines3= splice @lines2,0,27; #print "@lines2\n";

if ("normal_variant_filtered.vcf" eq "tumor_variant_filtered.vcf") { print "They're equal\n"; }else{ print "They are not equal\n"; }

foreach (@lines) { if ($_ =~ /^#/) { next; } else { my @stuff = split /\s/, $_; # print "$stuff[1]"; my $string = "1,$stuff[1],1,$stuff[3]/$stuff[4]\n"; $sift_string = $sift_string . $string; my @stuff2 = split /\s/, $_; # print "$stuff[1]"; #print "$sift_string\n"; my $string2 = "1,$stuff[1],1,$stuff[3]/$stuff[4]\n"; $sift_string2 = $sift_string . $string; print "$sift_string2\n"; } }

i tried doing something like this but didn't help.