text: user982: I have download a batch of refseq fasta files and want to split them based on strain. this is complicated by accession numbers, project IDs etc. I downloaded the genomes via eDirect from NCBI. I know I can split the output but I want them together for some preliminary parsing (remove ambiguous and plasmid sequences etc.)
  See example below. 
  >NZ_AESB01000237.1 Clostridium botulinum BKT028387 contig00237 ATGC
  >NZ_AESB01000236.1 Clostridium botulinum BKT028387 contig00236 ATGC
  >NZ_ALNK01000035.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290634 ATGC
  >NZ_ALNK01000031.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290635 ATGC
  
  I want two files:
  Clostridium botulinum BKT028387.fas and Peptoanaerobacter stomatis OBRC8.fas
  My actual dataset is much larger (~70,000 genomes) and not too sure about the headers in each case.
  
  Im thinking the following could work but Im not sure how to execute it:
  Using a grep command to extract the lines but split the output.
  If I download the summary list from refseq and use the names as the query example:
  grep -Fwf refseq.list sequences.fa
  
  
  Summary file
  
  Peptoanaerobacter stomatis strain OBRC8, whole genome shotgun sequencing project 7,810,226 rc other DNA  NZ_AJXR00000000.1 GI:938474759
  
  Clostridium botulinum BKT028387, whole genome shotgun sequencing project 4,134,593 rc other DNA  NZ_NOJY00000000.1 GI:1231817869
  
  
  But I could edit it down to 
  Peptoanaerobacter stomatis strain OBRC8
  Clostridium botulinum BKT028387
  
  Am I one the right track. How do I finish this?
  
  user2046: The SeqBuddy --pull_records command can grab sequences based on a regular expression against record IDs and/or metadata. Given your example:
  $: seqbuddy.py sequences.fa --pull_records "Clostridium botulinum.*BKT028387" full > Clostridium_botulinum_BKT028387.fas
  $: seqbuddy.py sequences.fa --pull_records "Peptoanaerobacter stomatis.*OBRC8" full > Peptoanaerobacter_stomatis_OBRC8.fas
  
  
  user2580: I had a similar problem and found a solution using awk. Since this is the first time I use awk, the syntax might not be perfect, but you should get an idea of what is possible.
  awk -F' ' '{if(/^>/){sub("strain ", ""); name=$2" "$3" "$4; print > name".fa"}else{print > name".fa"}}' yourfile.fa
  
  (more information below)
  This should transform a file like this:
  yourfile.fa
  >NZ_AESB01000237.1 Clostridium botulinum BKT028387 contig00237 ATGC
  ATGCATGC
  ATGCATGC
  >NZ_AESB01000236.1 Clostridium botulinum BKT028387 contig00236 ATGC
  ATGCATGC
  >NZ_ALNK01000035.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290634 ATGC
  ATGCATGC
  ATGCATGC
  >NZ_ALNK01000031.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290635 ATGC
  ATGCATGC
  
  Into two files like this:
  'Clostridium botulinum BKT028387.fa'
  >NZ_AESB01000237.1 Clostridium botulinum BKT028387 contig00237 ATGC
  ATGCATGC
  ATGCATGC
  >NZ_AESB01000236.1 Clostridium botulinum BKT028387 contig00236 ATGC
  ATGCATGC
  
  'Peptoanaerobacter stomatis OBRC8.fa'
  >NZ_ALNK01000035.1 Peptoanaerobacter stomatis OBRC8 ctg120007290634 ATGC
  ATGCATGC
  ATGCATGC
  >NZ_ALNK01000031.1 Peptoanaerobacter stomatis OBRC8 ctg120007290635 ATGC
  ATGCATGC
  
  How it works
  
  -F' ' tells awk to split the line by space (accessible with $1, $2, ...)
  if(/^>/) looks for line starting with the > character, then sub("strain ", "") remove the word strain and name=$2" "$3" "$4; print > name".fa" print the line into a new file with a name based on what you asked for.
  else if the line does not start with >, {print > name".fa"} appends it to a file using the current value of name.
  
  This solution will work regardless of the number of lines between the fasta headers.
  Note that I won't recommend using spaces in filenames, so I would replace name=$2" "$3" "$4 with something like name=$2"_"$3"_"$4.
  I used this solution to split a large fasta file (>5 GB) containing >5k genomes with >25k scaffolds and I found it fast and rather memory efficient.
  
original_text: 982: <p>I have download a batch of refseq fasta files and want to split them based on strain. this is complicated by accession numbers, project IDs etc. I downloaded the genomes <a href="https://bioinformatics.stackexchange.com/questions/3012/ncbi-edirect-suite-to-download-all-genome-sequences-associated-with-a-query">via eDirect from NCBI</a>. I know I can split the output but I want them together for some preliminary parsing (remove ambiguous and plasmid sequences etc.)</p>
  
  <p>See example below. </p>
  
  <pre><code>&gt;NZ_AESB01000237.1 Clostridium botulinum BKT028387 contig00237 ATGC
  &gt;NZ_AESB01000236.1 Clostridium botulinum BKT028387 contig00236 ATGC
  &gt;NZ_ALNK01000035.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290634 ATGC
  &gt;NZ_ALNK01000031.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290635 ATGC
  </code></pre>
  
  <p>I want two files:
  Clostridium botulinum BKT028387.fas and Peptoanaerobacter stomatis OBRC8.fas</p>
  
  <p>My actual dataset is much larger (~70,000 genomes) and not too sure about the headers in each case.</p>
  
  <hr>
  
  <p>Im thinking the following could work but Im not sure how to execute it:</p>
  
  <p>Using a <a href="https://unix.stackexchange.com/questions/361616/extract-lines-from-large-tab-delimited-file-using-a-list">grep command</a> to extract the lines but split the output.</p>
  
  <p>If I download the summary list from refseq and use the names as the query example:</p>
  
  <pre><code>grep -Fwf refseq.list sequences.fa
  </code></pre>
  
  <blockquote>
    <p>Summary file</p>
    
    <ol>
    <li><p>Peptoanaerobacter stomatis strain OBRC8, whole genome shotgun sequencing project 7,810,226 rc other DNA  NZ_AJXR00000000.1 GI:938474759</p>
    
    <ol start="2">
    <li>Clostridium botulinum BKT028387, whole genome shotgun sequencing project 4,134,593 rc other DNA  NZ_NOJY00000000.1 GI:1231817869</li>
    </ol></li>
    </ol>
    
    <p>But I could edit it down to </p>
    
    <p>Peptoanaerobacter stomatis strain OBRC8</p>
    
    <p>Clostridium botulinum BKT028387</p>
  </blockquote>
  
  <p>Am I one the right track. How do I finish this?</p>
   2046: <p>The <a href="https://github.com/biologyguy/BuddySuite/wiki/SB-Pull-records" rel="nofollow noreferrer">SeqBuddy --pull_records</a> command can grab sequences based on a regular expression against record IDs and/or metadata. Given your example:</p>
  
  <pre><code>$: seqbuddy.py sequences.fa --pull_records "Clostridium botulinum.*BKT028387" full &gt; Clostridium_botulinum_BKT028387.fas
  $: seqbuddy.py sequences.fa --pull_records "Peptoanaerobacter stomatis.*OBRC8" full &gt; Peptoanaerobacter_stomatis_OBRC8.fas
  </code></pre>
   2580: <p>I had a similar problem and found a solution using awk. Since this is the first time I use awk, the syntax might not be perfect, but you should get an idea of what is possible.</p>
  
  <pre><code>awk -F' ' '{if(/^&gt;/){sub("strain ", ""); name=$2" "$3" "$4; print &gt; name".fa"}else{print &gt; name".fa"}}' yourfile.fa
  </code></pre>
  
  <p><em>(more information below)</em></p>
  
  <p>This should transform a file like this:</p>
  
  <p>yourfile.fa</p>
  
  <pre><code>&gt;NZ_AESB01000237.1 Clostridium botulinum BKT028387 contig00237 ATGC
  ATGCATGC
  ATGCATGC
  &gt;NZ_AESB01000236.1 Clostridium botulinum BKT028387 contig00236 ATGC
  ATGCATGC
  &gt;NZ_ALNK01000035.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290634 ATGC
  ATGCATGC
  ATGCATGC
  &gt;NZ_ALNK01000031.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290635 ATGC
  ATGCATGC
  </code></pre>
  
  <p>Into two files like this:</p>
  
  <p>'Clostridium botulinum BKT028387.fa'</p>
  
  <pre><code>&gt;NZ_AESB01000237.1 Clostridium botulinum BKT028387 contig00237 ATGC
  ATGCATGC
  ATGCATGC
  &gt;NZ_AESB01000236.1 Clostridium botulinum BKT028387 contig00236 ATGC
  ATGCATGC
  </code></pre>
  
  <p>'Peptoanaerobacter stomatis OBRC8.fa'</p>
  
  <pre><code>&gt;NZ_ALNK01000035.1 Peptoanaerobacter stomatis OBRC8 ctg120007290634 ATGC
  ATGCATGC
  ATGCATGC
  &gt;NZ_ALNK01000031.1 Peptoanaerobacter stomatis OBRC8 ctg120007290635 ATGC
  ATGCATGC
  </code></pre>
  
  <p><strong>How it works</strong></p>
  
  <ul>
  <li><code>-F' '</code> tells awk to split the line by space (accessible with <code>$1</code>, <code>$2</code>, ...)</li>
  <li><code>if(/^&gt;/)</code> looks for line starting with the <code>&gt;</code> character, then <code>sub("strain ", "")</code> remove the word <code>strain</code> and <code>name=$2" "$3" "$4; print &gt; name".fa"</code> print the line into a new file with a name based on what you asked for.</li>
  <li>else if the line does not start with <code>&gt;</code>, <code>{print &gt; name".fa"}</code> appends it to a file using the current value of <code>name</code>.</li>
  </ul>
  
  <p>This solution will work regardless of the number of lines between the fasta headers.
  Note that I won't recommend using spaces in filenames, so I would replace <code>name=$2" "$3" "$4</code> with something like <code>name=$2"_"$3"_"$4</code>.
  I used this solution to split a large fasta file (>5 GB) containing >5k genomes with >25k scaffolds and I found it fast and rather memory efficient.</p>
  