0
$\begingroup$

I am trying to simulate metagenome sequencing. So I will start with a file with a lot of concatenated genomes. From there, I would like to randomly extract 10,000 sequences of length 200bps.

I got help subtracting a random 200 sequence from a fasta file like this:

>JQ086372.1 Enterobacteria phage HK446, complete genome 
ATGGCAGATAAACGAATCCGTTCCGACAGTTCGGCGGCAGCGGTTCAGGCCATGAAAAATGCAGCTGTGGCTTCTCCCGCTTCAGTCTGCCAATATGGATGTGAAACTCGTCGGAAAAAAAGTGGTTTATCGCTATCAACGCAACAGCGAATACGCCGACTTTTCGCAGAAAGAGATTTTTCACCTTAAAGGCTTCGGATTCACCGGGCTTGTAGGCCTGTCAATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTAC...(this still goes for a lot more A or C or G or T)

using this:

tail -n +2 file.fasta | tr -d '\n' > newfile
n=$(stat -c "%s" newfile)
r=$(shuf -i1-"$((n-200+1))" -n1)
newfile tail -c+"$r" | head -c200

And as an output I got: GTTAAGGGCGCCGTTCTGCTTGTCTTTGCTGACATGTTTGAACACCGGACGGCACAAAGCGAAGTGCAGCTTTATGAGAATGCAGCCGCCGAACGCATGATGTTCATTCATCGCAACTGGCGCGGTAAATCTGAACCTGAGGAGGGCTCCTGATGGAACCTGGACGATTCAGGCACAGGGTAAAAATTCTCACCTTCACG

But now I want to extract instead of 1 fragment of 200 characters; 10,0000 random fragments of 200 characters using the strategy above. For this, I am looking to exclude all those lines that begin with ">" (p.e. >JQ086372.1 Enterobacteria phage HK446, complete genome).The input file looks like this:

>JQ086372.1 Enterobacteria phage HK446, complete genome 
ATGGCAGATAAACGAATCCGTTCCGACAGTTCGGCGGCAGCGGTTCAGGCCATGAAAAATGCAGCTGTGGCTTCTCCCGCTTCAGTCTGCCAATATGGATGTGAAACTCGTCGGAAAAAAAGTGGTTTATCGCTATCAACGCAACAGCGAATACGCCGACTTTTCGCAGAAAGAGATTTTTCACCTTAAAGGCTTCGGATTCACCGGGCTTGTAGGCCTGTCAATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTAC...(this still goes for a lot more A or C or G or T)

>NC_010391.1 Salmonella phage Fels-1, complete genome 
TGATTGCATAACCACTTAACATCTTGTTTTATCTAAATAAAATTAAGCATGTTATCTTTTTGGGGCACTCCTGGGGCAGTAGATGCCAGTTGTTGATTCAGTATATCTACTTGTGACTGGTTATTATCCCGATTTTTTTAGTTTTAAGGTGTTGACATAGCCATCCATGCTCCATATACTGTATAGACCATCTGAGCGTTTGAGTGCCCCATTTGTGAAGCAATAAAGTTCGGGTTCGCGCCAGCGGCAAGCGCCCAGCATGCACCGATTTTTTTAGTTTTAAGGTGTTGACATTAGGTATGTCGGGACTGGTATGCTTTCCTGTGTCGCAGCCCGGCGCGTCTCAATGCAGATTCCCATATCCTGTTCAGGGAAGGTGCAGCATAATGTGCTTTGGATTTACCTGAGCGATCGGTAATTTGCGGATTGAAGACAAAGGTAATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTAC...(this still goes for a lot more A or C or G or T)

>NC_007577.1 Prochlorococcus marinus str. MIT 9312, complete genome 
AATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTACAATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTAC

So for my new output, I will end up with a file with 10,000 subsets of 200 characters that would look like this:

>1
ATGGCAGATAAACGAATCCGTTCCGACAGTTCGGCGGCAGCGGTTCAGGCCATGAAAAATGCAGCTGTGGCTTCTCCCGCTTCAGTCTGCCAATATGGATGTGAAACTCGTCGGAAAAAAAGTGGTTTATCGCTATCAACGCAACAGCGAATACGCCGACTTTTTCAGTCTGCCAATATGCGTTCCGACAGCAGATAAAC

>2
TGATTGCATAACCACTTAACATCTTGTTTTATCTAAATAAAATTAAGCATGTTATCTTTTTGGGGCACTCCTGGGGCAGTAGATGCCAGTTGTTGATTCAGTATATCTACTTGTGACTGGTTATTATCCCGATTTTTTTAGTTTTAAGGTGTTGACATAGCCATCCATGCTCCATATACTGTATAGACCATCTGCAGATA

>3
AATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAAT

Makes sense? Any help?

Thanks in advance

$\endgroup$
0
0
$\begingroup$

Given that you want to simulate short-read sequencing from reference genomes, there are read simulator programs that already exist. For example, see this answer.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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