Thank you, I'll try to do something with it.
Dear colleagues, I am applying for your help once more. Here is my perl script for counting the length between the binding sites and the start points of exons. So, there are some genes saved in the file called sequence.txt and some amount of binding sites in motif.txt. The thing I want to do is to count the length of the fragment for every gene if it has any of the binding sites from motif.txt. This current script does not work in the way I would like it. How should it be changed? Shall I add another while loop for $motif?
$string_filename = 'sequence.txt';
open(FILE, $string_filename) || die("Couldn't read file $string_filename\n");
$motif_filename = 'motif.txt';
open(MOTIF, $motif_filename) || die("Couldn't read file $motif_filename\n");
local $/ = "\n>";
while (my $seq = <FILE>) {
chomp $seq;
$seq =~ s/^>*.+\n//;
$seq =~ s/\n//g;
$R = length $seq;
$motif = <MOTIF>;
chomp $motif;
$motif =~ s/^>*.+\n//;
$motif =~ s/\n//g;
if ( $seq =~ /$motif/ ) { ## insert actual binding site
$M = $';
$W = length $M;
if ( $seq =~ /[A-Z]/) { ## exon start
$K = $`;
$Z = length $K;
$x = $W + $Z - $R;
print "\n\ the distance is the following: $x\n\n";
} else {
print "\n\ I couldn't find the start codon.\n\n";
}
} else {
print "\n\ I couldn't find the binding site.\n\n";
}
}
close MOTIF;
close FILE;
exit;
4 answers
I think you need to take a step back and think about the logic of your problem, before writing the code. Break your problem down into steps, like this:
- Read motifs into a data structure over which you can loop (e.g. an array)
- Read in your sequence(s)
- For each sequence, identify segments between binding site and exon start
- For each segment...
- For each motif, check for match to segment
- If at least one match, store/print segment length (and perhaps other useful information)
For tools to convert FASTQ-FASTA, see comments under your question.
As we said in your previous query, Bioperl makes this kind of task very trivial. There is a learning curve, but I would urge you to investigate - it will pay off in the long term. Using more "standard" code also makes it a lot easier for other people to read and debug.
you should load your motifs into an array and then iterate over that motif array for each sequence
I eventually recognised how to modify my script to do the job I wished. Now it looks like the following:)
$string_filename = 'seq.txt';
open(FILE, $string_filename) || die("Couldn't read file $string_filename\n");
$motif_filename = 'motif.txt';
open(MOTIF, $motif_filename) || die("Couldn't read file $motif_filename\n");
local $/ = "\n>";
while (@motif = <MOTIF>) {
while (my $seq = <FILE>) {
chomp $seq;
$seq =~ s/^>*.+\n//;
$seq =~ s/\n//g;
$R = length $seq;
foreach $site (@motif) {
chomp $site;
$site =~ s/^>*.+\n//;
$site =~ s/\n//g;
if ( $seq =~ /$site/ ) {
$M = $';
$W = length $M;
if ( $seq =~ /[AGTC]/) {
$K = $`;
$Z = length $K;
$x = $W + $Z - $R;
print "\nthe distance is the following: $x\n\n";
}
}
}
}
}
close FILE;
exit;
Today I tried to do the same using BioPerl module and recognised it to be much more easier indeed. Thanks a million for your help and support:)
accttttgat gaagatcagc 2221 atacacaaat tacaaaagtc tgaatttttt tttatcaaga gggataaaac accatgaaaa 2281 taaacttgaa taaactgaaa atggaccttt ttttttttaa tggcaatagg acattgtgtc 2341 agattaccag ttataggaac aattctcttt tcctgaccaa tcttgtttta ccctatacat 2401 ccacagggtt ttgacacttg ttgtccagtt gaaaaaaggt tgtgtagctg tgtcatgtat 2461 ataccttttt gtgtcaaaag gacatttaaa attcaattag gattaataaa gatggcactt 2521 tcccgtttta ttccagtttt ataaaaagtg gagacagact gatgtgtata cgtaggaatt 2581 ttttcctttt gtgttctgtc accaactgaa gtggctaaag agctttgtga tatactggtt 2641 cacatcctac ccctttgcac ttgtggcaac agataagttt gcagttggct aagagaggtt 2701 tccgaagggt tttgctacat tctaatgcat gtattcgggt taggggaatg gagggaatgc 2761 tcagaaagga aataatttta tgctggactc tggaccatat accatctcca gctatttaca 2821 cacacctttc tttagcatgc tacagttatt aatctggaca ttcgaggaat tggccgctgt 2881 cactgcttgt tgtttgcgca ttttttttta aagcatattg gtgctagaaa aggcagctaa 2941 aggaagtgaa tctgtattgg ggtacaggaa tgaaccttct gcaacatctt aagatccaca 3001 aatgaaggga tataaaaata atgtcatagg taagaaacac agcaacaatg acttaaccat 3061 ataaatgtgg aggctatcaa caaagaatgg gcttgaaaca ttataaaaat tgacaatgat 3121 ttattaaata tgttttctca attgtaacga cttctccatc tcctgtgtaa tcaaggccag 3181 tgctaaaatt cagatgctgt tagtacctac atcagtcaac aacttacact tattttacta 3241 gttttcaatc ataatacctg ctgtggatgc ttcatgtgct gcctgcaagc ttcttttttc 3301 tcattaaata taaaatattt tgtaatgctg cacagaaatt ttcaatttga gattctacag 3361 taagcgtttt ttttctttga agatttatga tgcacttatt caatagctgt cagccgttcc 3421 acccttttga ccttacacat tctattacaa tgaattttgc agttttgcac attttttaaa 3481 tgtcattaac tgttagggaa ttttacttga atactgaata catataatgt ttatattaaa 3541 aaggacattt gtgttaaaaa ggaaattaga gttgcagtaa actttcaatg ctgcacacaa 3601 aaaaaagaca tttgattttt cagtagaaat tgtcctacat gtgctttatt gatttgctat 3661 tgaaagaata gggttttttt tttttttttt tttttttttt ttaaatgtgc agtgttgaat 3721 catttcttca tagtgctccc ccgagttggg actagggctt caatttcact tcttaaaaaa 3781 aatcatcata tatttgatat gcccagactg catacgattt taagcggagt acaactacta 3841 ttgtaaagct aatgtgaaga tattattaaa aaggtttttt tttccagaaa tttggtgtct 3901 tcaaattata ccttcacctt gacatttgaa tatccagcca ttttgtttct taatggtata 3961 aaattccatt ttcaataact tattggtgct gaaattgttc actagctgtg gtctgaccta 4021 gttaatttac aaatacagat tgaataggac ctactagagc agcatttata gagtttgatg 4081 gcaaatagat taggcagaac ttcatctaaa atattcttag taaataatgt tgacacgttt 4141 tccatacctt gtcagtttca ttcaacaatt tttaaatttt taacaaagct cttaggattt 4201 acacatttat atttaaacat tgatatatag agtattgatt gattgctcat aagttaaatt 4261 ggtaaagtta gagacaacta ttctaacacc tcaccattga aatttatatg ccaccttgtc 4321 tttcataaaa gctgaaaatt gttacctaaa atgaaaatca acttcatgtt ttgaagatag 4381 ttataaatat tgttctttgt tacaatttcg ggcaccgcat attaaaacgt aactttattg 4441 ttccaatatg taacatggag ggccaggtca taaataatga cattataatg ggcttttgca 4501 ctgttattat ttttcctttg gaatgtgaag gtctgaatga gggttttgat tttgaatgtt 4561 tcaatgtttt tgagaagcct tgcttacatt ttatggtgta gtcattggaa atggaaaaat 4621 ggcattatat atattatata tataaatata tattatacat actctcctta ctttatttca 4681 gttaccatcc ccatagaatt tgacaagaat tgctatgact gaaaggtttt cgagtcctaa 4741 ttaaaacttt atttatggca gtattcataa ttagcctgaa atgcattctg taggtaatct 4801 ctgagtttct ggaatatttt cttagacttt ttggatgtgc agcagcttac atgtctgaag 4861 ttacttgaag gcatcacttt taagaaagct tacagttggg ccctgtacca tcccaagtcc 4921 tttgtagctc ctcttgaaca tgtttgccat acttttaaaa gggtagttga ataaatagca 4981 tcaccattct ttgctgtggc acaggttata aacttaagtg gagtttaccg gcagcatcaa 5041 atgtttcagc tttaaaaaat aaaagtaggg tacaagttta atgtttagtt ctagaaattt 5101 tgtgcaatat gttcataacg atggctgtgg ttgccacaaa gtgcctcgtt tacctttaaa 5161 tactgttaat gtgtcatgca tgcagatgga aggggtggaa ctgtgcacta aagtgggggc 5221 tttaactgta gtatttggca gagttgcctt ctacctgcca gttcaaaagt tcaacctgtt 5281 ttcatataga atatatatac taaaaaattt cagtctgtta aacagcctta ctctgattca 5341 gcctcttcag atactcttgt gctgtgcagc agtggctctg tgtgtaaatg ctatgcactg 5401 aggatacaca aaaataccaa tatgatgtgt acaggataat gcctcatccc aatcagatgt 5461 ccatttgtta ttgtgtttgt taacaaccct ttatctctta gtgttataaa ctccacttaa 5521 aactgattaa agtctcattc ttgtcaaaaa aaaaaaaaaa aaaaaaaaaa aa
Log in to answer this question.
please, show us what your inputs look like.
the binding sites look like this
and the genes
binding site=A FASTQ file ??!!!!!
I recognise that it is a problem but i haven't found the way to convert it to fasts yet.
There is a FASTQ-FASTA converter here - and Bioperl's
Bio::SeqIO::fastqwill handle some FASTQ formats.And further discussion of FASTQ-FASTA conversion at StackOverflow.