I am huge supporter of Python as a language for beginner over PERL for various reasons:
https://www.linkedin.com/pulse/perl-verus-python-high-time-first-language-beginners-vijay-lakhujani
I want to match my motif file to the set of sequences and get the desired output like this:
AAAAAA 2
AAAAAT 1
AAAAAG 1
AAAAAC 0
the perl code I'm using is
open(FILE,"dataset1.txt");
%jobs_list = <FILE>;
@indexes=keys%jobs_list;
@seq=values%jobs_list;
open(IN,"Motif6.txt");
@cis = <IN>;
$count=0;
for ($i=0;$i<=$#seq;$i++) {
foreach $index (@cis) {
if($seq[$i] =~/$index/) {
$count++;
}
}
print "$index\t$count \n";
}
the data set file is -
>1
AAAAAAATGCATCGATCGAC
>2
AAAAAAGTCGATCGATCAGC
and the Motif6 file is-
AAAAAA
AAAAAT
AAAAAG
AAAAAC
Please suggest me what can I do to make it correct and get the desired output.
Perl is difficult to learn for beginners without any programming experiment.
Maybe we can paste some code of some basic tasks like reading FASTA sequences to give some examples and encourage beginners to make the first step.
Reading code of others is also an important way of learning programming language.
It's really difficult to debug others' Perl code. Here's a rewritten version with detailed comments.
Seqs:
$ cat dataset1.txt
>1
AAAAAAATGC
ATCGATCGAC
>2
AAAAAAGTCG
ATCGATCAGC
Motifs:
$ cat Motif6.txt
AAAAAA
AAAAAT
AAAAAG
AAAAAC
Code:
#!/usr/bin/env perl
use strict;
# read motifs
my $file_motif = "Motif6.txt";
my @motifs = (); # an array to store motifs
open my $fh, "<", $file_motif or die "fail to open file: %file_motif";
for my $line (<$fh>) {
$line =~ s/\r?\n//g; # chomp "\r" and "\n"
next if $line eq "" or $line =~ m/\s+/; # skip blank line
push @motifs, $line; # append to array
}
close $fh;
# read sequences and count motifs
my $file_seq = "dataset1.txt";
my %counter = (); # a hashmap to store countings
open my $fh, "<", $file_seq or die "fail to open file: $file_seq";
my ($header, $seq) = ("", "");
for my $line (<$fh>) {
$line =~ s/\r?\n//g; # chomp "\r" and "\n"
if (substr( $line, 0, 1 ) eq '>') { # FASTA header line
unless ($header eq "" and $seq eq "") { # previous sequence
&count_motifs($header, $seq);
}
$header = substr ($line, 1); # update header
$seq = ""; # reset seq
} else { # sequence line
$seq = $seq . $line # concatenate sequence
}
}
close $fh;
unless ($header eq "" and $seq eq "") { # do not forget the last sequence
&count_motifs($header, $seq);
}
# print result, in decending order of count
for my $motif (sort {$counter{$b} <=> $counter{$a}} keys %counter) {
print "$motif\t$counter{$motif}\n";
}
sub count_motifs {
my ($header, $seq) = @_;
for my $motif (@motifs) {
my $n = 0;
my $len_motif = length $motif;
my ($begin, $end) = (-1, -1);
# if you just want to test whether the seq contains a motif, use:
$n++ if $seq =~ m/$motif/i;
# if you want to find all matches in sequences, use:
# while ($seq =~ m/$motif/gi) { # use regular expression to locate motifs
# my $pos = pos $seq; # http://perldoc.perl.org/functions/pos.html
# ($begin, $end) = ($pos - $len_motif + 1, $pos);
# print STDOUT "seq: $header\tmotif: $motif\tlocation: $begin-$end\n";
# pos $seq = $pos - $len_motif + 1;
# $n++;
# }
$counter{$motif} += $n;
}
}
Result:
$ perl count_motifs_in_seqs.pl
AAAAAA 2
AAAAAG 1
AAAAAT 1
AAAAAC 0
I am huge supporter of Python as a language for beginner over PERL for various reasons:
https://www.linkedin.com/pulse/perl-verus-python-high-time-first-language-beginners-vijay-lakhujani
So do I, i alway advise self-learning beginner choose python not perl. i use python most of the time after quiting perl.
Some comments:
use strict; use warningschomp, and to skip blank lines. That also applies to your Motif file where you have blank lines, and they will be searched in the sequence.foreach or map Logic of counting and iteration does not combine, you defined only a global counter, but your loop logic does not fit to this. You can easily see that the following code will output an ever increasing number for each combination of motif and target sequence:
# this is pseudocode of what you are trying:
0: counter := 0
1: forall target in target sequences
2: forall motif in motifs
if (motif matches target):
increment (counter)
print "motif occurs in" counter "sequences at least once (NOT)"
You need to change the sequence of loop control in line 0-2 to suit your problem (hint: 2,0,1 looks like a good sequence)
Thankyou sir for your suggestion it helped me a big tym!
Check this post counting motifs in dna sequence, no programming needed. You might skip the step 1.
I am guessing OP just wants or needs to learn perl, or it is some kind of homework.
ok, let's leave him/her alone :)
I think if someone know the answer then he should tell that even if he knows that someone is asking for homework, it can be beneficial in many ways
I disagree. We're all volunteers here and we are never obliged to give an answer. Furthermore, homework is supposed to learn you something. In my opinion, the best way to learn how to program is to struggle and fight with the code until it does what you want to. Building confidence to use these sites by answering homework questions seems like a pointless way, besides, we don't need to market. If you like biostars, you are welcome. But you are not a customer.
Your last reason might be true occasionally, but definitely not for this question. This text and title are far too aspecific.
I also disagree. In addition to providing practical experience in coding, homework also provides the instructor feedback on your mastery of the material. This feedback allows him/her to provide additional assistance when needed. It also allows the assessment of a grade/score for the course, which is an (admittedly rough) indicator of your proficiency to third parties.
You can use hash to store your motif,and the values are the numbers of motif.If match one motif,Its value should +1.
#! /usr/bin/env perl
open IN,"PATH" || die "$!";
#AAAAAAATGCATCGATCGAC
#AAAAAAGTCGATCGATCAGC
open MOTIF,"PATH" || die "$!";
#AAAAAA
#AAAAAT
#AAAAAG
#AAAAAC
%motif=();
while(defined($t=<MOTIF>)){
chomp($t);
$motif{$t}=0;
}
while(defined($line=<IN>)){
chomp($line);
foreach (keys %motif){
my $count=$line=~m/$_/;
$motif{$_}+=$count;
}
}
@motif=keys %motif;
foreach (@motif){
print "$_\t$motif{$_}\n";
}
close IN;
close MOTIF;
And the result:
AAAAAG 1
AAAAAT 1
AAAAAA 2
AAAAAC 0
Log in to answer this question.
I've done what I can to fix the formatting.
how i can see the formatting?
Start with this:
Edit: And don't forget that when reading lines from a file, you're also reading the end of line character(s).
Do you want to do it in PERL only? you can get this done using simple
grepcommandsir i have to do it in perl only!!!
i see, it's a homework:)
An entire separate discussion, but is it still 'valid' to 'force' students to use perl? Okay the language was hugely important during the human genome project, but aren't those times changing?
I am not sure we should start with this discussion again. Perl is an extremely stable, high performance and well debugged programming language with the most comprehensive set of Bio::* libraries. Large pipelines such as miRdeep and Trinity are implemented in Perl.
Absolutely. For that change to go to completion younger folk like you need to get their PhD's quickly and replace old perl users.
On a more serious note if you are learning on your own then you are free to choose the latest and greatest but if you are taking a class then you don't have much choice.
Fast cheap hardware/memory has made programmer's job easy. See what Margaret Hamilton had to work with.
sir its not a homework its a part of some research work in which i was stuck upon.
Then why did you absolutely have to do it in perl?
Maybe Perl is his main programming languages.
Absolutely,bcoz the thing m working on I started on with perl and dis was a part of it I had to continue in perl programming only! Anyways thankyou all for helping through it!
You should have made your question more clear, so you an get help more quickly. If you tell this at the very beginning, I can tell you this simplest way, which can save you and others' time:
Using my lovely SeqKit and csvtk:
Step 1. convert the list file containing motif sequences to FASTA format
Step 2. locate motifs and count
You can do this in Windows/Linux/Mac, because both seqkit and csvtk support them.
Oh m sorry next time if I get stuck ill try to make questions more clear!:)