Extract match and everything before next match, do it for every match

Question

I have a file with this structure:

>Cluster 0
0       51aa, >MG00HS05:520:C8M1TACXX:3:1101:1428:2080/1... *
1       51aa, >MG00HS05:520:C8M1TACXX:3:1101:1658:2480/1... at 3:51:1:49/96.08%
2       51aa, >MG00HS05:520:C8M1TACXX:3:1101:15131:2756/1... at 1:51:1:51/100.00%
[thousands of similarly looking lines]
>Cluster 1
0       51aa, >MG00HS05:520:C8M1TACXX:3:1101:3733:2088/1... *
1       50aa, >MG00HS05:520:C8M1TACXX:3:1101:6962:2026/1... at 2:50:1:49/98.00%
2       51aa, >MG00HS05:520:C8M1TACXX:3:1101:14617:2071/1... at 2:51:1:50/96.08%
[thousands of similarly looking lines]
>Cluster 2
0       51aa, >MG00HS05:520:C8M1TACXX:3:1101:5164:2153/1... *
1       51aa, >MG00HS05:520:C8M1TACXX:3:1101:15660:20057/1... at 1:51:1:51/98.04%
2       51aa, >MG00HS05:520:C8M1TACXX:3:1101:8563:35493/1... at 1:50:1:51/96.08%
[thousands of similarly looking lines]

The lines starting with > are about two millions.

I'd like to extract the lines starting with > and those following it, without taking the following line starting with >, and put them in files. Something like this:

File_one:

>Cluster 0
0       51aa, >MG00HS05:520:C8M1TACXX:3:1101:1428:2080/1... *
1       51aa, >MG00HS05:520:C8M1TACXX:3:1101:1658:2480/1... at 3:51:1:49/96.08%
2       51aa, >MG00HS05:520:C8M1TACXX:3:1101:15131:2756/1... at 1:51:1:51/100.00%
[thousands of similarly looking lines]

File_two

>Cluster 1
0       51aa, >MG00HS05:520:C8M1TACXX:3:1101:3733:2088/1... *
1       50aa, >MG00HS05:520:C8M1TACXX:3:1101:6962:2026/1... at 2:50:1:49/98.00%
2       51aa, >MG00HS05:520:C8M1TACXX:3:1101:14617:2071/1... at 2:51:1:50/96.08%
[thousands of similarly looking lines]

File_three

>Cluster 2
0       51aa, >MG00HS05:520:C8M1TACXX:3:1101:5164:2153/1... *
1       51aa, >MG00HS05:520:C8M1TACXX:3:1101:15660:20057/1... at 1:51:1:51/98.04%
2       51aa, >MG00HS05:520:C8M1TACXX:3:1101:8563:35493/1... at 1:50:1:51/96.08%
[thousands of similarly looking lines]

I've written a script that should do it in bash, but it's not working. I'm not a pro in bash scripting.

mkdir FemaleMito1_clusters
while read i
        do $i > FemaleMito1_clusters/FemaleMito1_${i#>}
        n=1
        while [ `grep -A $n $i FemaleMito1_cdhit2 | tail -n1 | grep -c "^>"` -eq 0 ]
                do grep -A"$n" $i FemaleMito1_cdhit2 | tail -n1 >> FemaleMito1_clusters/FemaleMito1_"${i#>}"
                ((n++))
                done
        done < FemaleMito1_cdhit2_list #this is a file containing just the lines starting with >

How can I do it? Feel free to skip completely my script, there's probably a one-liner that does what I'm looking for.

I also have to filter the files and retain only the ones with more than a certain line number. I thought about doing it with a simple wc -l after creating the files, but if there's a way to include this in the command without creating useless files it's better.

Answer 1

You can do this quite easily in awk:

awk '{ if(/^>/){name=$0; sub(/^>/,"", name);}{print >> name".fa"}}' file.fa 

That will iterate over all of the input file's line and, if the first character is a >, it will save that line as name. Then, it will remove the > from the contents of name since you don't want that in the file name. Finally, each line is appended to a file called name.fa where name is whatever the current sequence's name is.

If you only want to print those sequences with more than N lines, you can use:

awk -v min=4 '{ 
               if(/^>/){ 
                    if(num >= min){
                        print seq >> name".fa"
                    } 
                    name=$0; 
                    sub(/^>/,"", name); 
                    seq=$0; 
                    num=0
                }
                else{
                    seq = seq"\n"$0; 
                    num++
                }
               }
               END{
                 if(num >= min){
                    print seq >> name".fa"
                 }
               }' file.fa 

As a general rule, don't use shell loops for text processing. They're slow, cumbersome, and prone to error.

Answer 2

Although (as you have been advised in comments) there are probably bioinformatics tools that may be more suitable for your application, it can be done using csplit:

csplit -sz file '/^>/' '{*}'

gives

$ head xx*
==> xx00 <==
>Number_one
[some thousands lines]

==> xx01 <==
>Number_two
[some other thousands lines, less than the latter]

==> xx02 <==
>Number_three
[Some other hundreds lines]

For options concerning the numbering and format of the output file names, refer to the manual page (man csplit)