Remove N's and split contigs

[Dagga]

Hi,

I have some genomes that I will be uploading to NCBI soon. I have been told that all N's need to be removed and the contigs split at this position.

I am new to command line interface so I was hoping someone could recommend a program and simple script that could do this for me. I would like to remove all N's and then split the contig at the location of the N's results in two new contigs. For example

Contig 1: ATCGGATAANNNNNNNNNATCGCCGAT

Contig 1.1: ATCGGATAA

Contig 1.2 ATCGCCGAT


Thanks!

[TiborNagy]

perl -ne 'if($_ =~ /([^N]+)N+([^N]+)/){print $1;print stderr $1}' input.seq >contig1.txt 2>contig2.txt

It will split the input file (input.seq) into contig1.txt and contig2.txt

[mastal]

should that be

print stderr $2

[Dagga]

Thanks for that!!

Will this rename the contigs?

Will the contig that is split be called the same thing in contig1.txt and contig2.txt.

Is it possible to rename the contigs when they are split. For example, if contig 84 is split into two contigs can they be renamed contig 84.1 and contig 84.2 for each half, respectively?

[TiborNagy]

mastal: you are right!
Dagga: This script does not handle the contig names, only sequences, because you do not tell us what kind of input format do you have.

[Dagga]

TiborNagy: Sorry, the file will be in fasta format post de novo assembly.

would you be able to alter the script to handle contig names please?

Thanks!

[mastal]

If you are doing your assemblies with velvet, setting '-scaffolding no' will stop velvet from joining contigs together with stretches of Ns.

[Dagga]

Excellent!

Whilst this does help with some genomes that I am assembling right now, we have some older genomes that were sequenced by BGI and these contain N's that we still need to have removed...

[TiborNagy]

Just for you :-)

#!/usr/bin/perl

$seq = "";

while(<>){
   chomp;

   if(/^>/){
      if($seq ne ""){
         if($seq =~ /([^N]+)N+([^N]+)/){
            print  "$id.1\n$1\n";
            print STDERR "$id.2\n$2\n";
         }
      }
      $seq = "";
      $id = $_;
   }
   else{
      $seq .= $_;
   }
}

if($seq =~ /([^N]+)N+([^N]+)/){
  print "$id.1\n$1\n";
  print STDERR "$id.2\n$2\n";
}

[Dagga]

Thanks!! I appreciate it!