1
$\begingroup$

I have seven fasta files, one per gene, with more than 400 fasta entries per file, like this: Input:

Gene1.fasta

>1721_1 gene name  
ATG   
>10_14 gene name  
GTT

Gene2.fasta

>1721_1 gene name  
TGA  
>10_14 gene name   
GAT

Output:

>1721_1  
ATGTGA  
>10_14   
GTTGAT

What I want to do is create one fasta (continuously) file containing seven genes per sample. As you can see, the sample name is not the same before the _ in the fasta headers.

Is there any way to do this with awk or other command-line tools? Thank you very much.

$\endgroup$
4
  • 1
    $\begingroup$ cat Gene*.fasta > all_genes.fasta ? $\endgroup$ Commented Oct 4, 2021 at 3:46
  • $\begingroup$ @TimurShtatland Thank you. But it should contain all genes sequences of each sample. Each sample should be separate files. $\endgroup$
    – Saruul B
    Commented Oct 4, 2021 at 11:25
  • $\begingroup$ I would suggest using samtools faidx to pull the sample's gene out of each file, and then removing the header. It will ignore anything after whitespace so gene name after space should not be problem. You could iterate across sample names in a bash loop. Of course, it requires that every sample have every gene. $\endgroup$ Commented Oct 5, 2021 at 11:57
  • $\begingroup$ okay. I will try $\endgroup$
    – Saruul B
    Commented Oct 11, 2021 at 4:31

3 Answers 3

2
$\begingroup$

The concat command from seqkit could concatenate sequences with same ID from multiple files.

seqkit concat [file1] [file2] ...
$\endgroup$
0
$\begingroup$

I think it can be done using python dictionary. You can split the file using the header and then make the fasta id as your key and the sequence as your value. Then using simple loop and conditions you can concatenate them.

$\endgroup$
2
  • $\begingroup$ Okay. I will try $\endgroup$
    – Saruul B
    Commented Oct 11, 2021 at 4:31
  • 1
    $\begingroup$ Hi and welcome to the site! Please try and post complete answers. This isn't really answering, it is only giving a general direction in which an answer could be found. $\endgroup$
    – terdon
    Commented Nov 5, 2021 at 17:47
0
$\begingroup$

Assuming your files are small enough to allow you to fit one of them into memory, you can do this:

$ awk '{ 
        if(/^>/){ name=$0 }
        else{ 
            if(NR==FNR){seq[name]=$0 } 
            else{
                print name; 
                printf "%s%s\n", seq[name],$0
            }
        }
      }' file1.fa file2.fa 
>1721_1 gene name
TGA  
>10_14 gene name
GTTGAT

If you cannot fit the files in memory, you can try this:

$ join -o1.1,1.2,2.2 -t $'\t' <(FastaToTbl file1.fa ) <(FastaToTbl file2.fa) | sed 's/\t//2' | TblToFasta 
>1721_1 gene name 
ATGTGA
>10_14 gene name 
GTTGAT

Note that this approach requires the FastaToTbl and TblToFasta scripts I have posted here: https://bioinformatics.stackexchange.com/a/2651/298.

Or, if brevity is more your thing, the same approach but more compact:

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.