How to get length of selected contigs within a fasta file generated by transcriptome assembly?

Question

[seta]

Dear all,

I would like to have length of some selected contigs within a fasta file resulted from de novo transcriptome assembly. I have names of selected contigs in the txt file like here, could you please help me out to get it?

>contig1
>contig2
>contig3

Thanks

[Rohit]

One of the easiest ways for this is to use samtools faidx contigs.fasta. The .fai file has the necessary information

If you need to extract prior then, remove the ">" and -

xargs samtools faidx ref.fasta < list >contigs.fasta

Answer 1

[Sej Modha]

Calculate length of each contig using following script and save the output of this script to an output file.

#!/usr/bin/perl -w

use strict;
#contig_length_fasta.pl

my $largest = 0;
my $contig = '';

if (@ARGV != 1) {
    print "contig_length_fastq fasta\n\n";
    exit;
}

my $filenameA = $ARGV[0];

open (IN, "$filenameA") or die "oops!\n";

my $read_name = '';
my $read_seq = '';

while (<IN>) {
    if (/^>(\S+)/) {
        $read_name = $1;
        $read_seq = "";

        while (<IN>) {
            if (/^>(\S+)/) {
                print "$read_name\t" . length($read_seq) . "\n";
                $read_name = $1;
                $read_seq = "";
            }
            else {
                chomp;
                $read_seq .= $_;
            }
        }
    }
}

close(IN);

print "$read_name\t" . length($read_seq) . "\n";

contig_length_fasta.pl input_file.fa > contig_length_out

Use simple grep to get length of the contigs of interest.

grep -f my_contig_list contig_length_out

Answer 2

[Matt Shirley]

Use pyfaidx:

$ (sudo) pip install pyfaidx
$ xargs faidx --transform chromsizes contigs.fasta < contigs.txt

Answer 3

[Lesley Sitter]

Hi,

I just adapted a short python script I have for these types of things. Is this something you want? It produces a file with the contigs, a tab and then the length of the sequences.

It takes two user arguments, the first is the fasta file, the second the name of a output file.

So just start it like this on a linux command line (just store the script as script.py or something)

$ python script.py fasta.fa seq_length.txt

#!usr/bin/env python

import sys 

file_open = open(sys.argv[1],'r').readlines() 
file_out = open(sys.argv[2],'w') 

switch = 0 
for lines in file_open:
    if '>' in lines:
        if switch == 0:
            file_out.write(lines.strip() + "\t")
            seq_length = 0
            switch = 1
        elif switch == 1:
            file_out.write(str(seq_length) + "\n" + lines.strip() + "\t")
            seq_length = 0
    else:
        seq_length+=len(lines.strip())
file_out.write(str(seq_length))

[seta]

Many thanks for all your comments. I try them

Answer 4

[andrew.j.skelton73]

Have you looked at this?

Answer 5

[5heikki]

If there are no linebreaks in sequences then, e.g.

awk '/^>contig1/{getline; print length($0)}' file.fasta