Commit d8f373ee authored by Anna Buschart's avatar Anna Buschart
Browse files

bug in getValidGenes1.R

parent 1b584a70
#!/usr/bin/perl
# this script takes a fasta file and a space-delimited file with IDs present in the fasta headers and start- and end- coordinates (with the first letter counted as position 1)
# 2 inputs: - the fasta file (eg. mOTU.v1.padded)
# - a space-separated table with IDs and coordinates (eg. mOTU.v1.padded.coord)
# 1 output to stdout: a fasta file with sequences indicated by the table
# Anna Heintz-Buschart, February 2015
use strict;
use warnings;
use Bio::DB::Fasta;
my $fastaFile = shift;
my $queryFile = shift;
my $db = Bio::DB::Fasta->new( $fastaFile );
open (IN, $queryFile);
while (<IN>){
chomp;
my @fields = split /\t/;
my $seq = $fields[0];
my $sequence = $db->seq($seq);
if (!defined( $sequence )) {
print STDERR "Sequence $seq not found. \n";
next;
}
my $start=0;
my $stop=0;
if ($fields[1] < $fields[2]) {
$sequence = substr($sequence,$fields[1]-1,($fields[2]-$fields[1]+1));
$start = $fields[1];
$stop = $fields[2];
} else {
$sequence = substr($sequence,$fields[2]-1,($fields[1]-$fields[2]+1));
$start = $fields[2];
$stop = $fields[1];
$sequence = join('', reverse(split(//, $sequence)));
$sequence =~ tr/ACGT/TGCA/;
}
print ">$seq\n", "$sequence\n";
}
close(IN);
...@@ -19,7 +19,7 @@ phyl2 <- phyl[ phyl$mOTU.species.annotation %in% pres$V1,] ...@@ -19,7 +19,7 @@ phyl2 <- phyl[ phyl$mOTU.species.annotation %in% pres$V1,]
mp <- merge(map,phyl2,by.x="mOTU",by.y=1) mp <- merge(map,phyl2,by.x="mOTU",by.y=1)
coord <- read.delim(paste("../../",COG,"/",COG,".mOTU.v1.padded.coord",sep=""),header=F,stringsAsFactors=F) coord <- read.delim(paste("../../",COG,"/",COG,".mOTU.v1.padded.coord",sep=""),header=F,stringsAsFactors=F)
co <- gsub(" ",".",coord$V1) co <- gsub(" ",".",coord$V1)
coord <- read.delim(paste("../../",COG,"/",COG,".mOTU.v1.padded.coord",sep=""),header=F,stringsAsFactors=F,sep=" ") #coord <- read.delim(paste("../../",COG,"/",COG,".mOTU.v1.padded.coord",sep=""),header=F,stringsAsFactors=F,sep=" ")
coord$name <- co coord$name <- co
mp2 <- merge(mp[mp$MarkerGene ==COG,],coord[,c(1,4)],by.x="SequenceID",by.y=1) mp2 <- merge(mp[mp$MarkerGene ==COG,],coord[,c(1,4)],by.x="SequenceID",by.y=1)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment