[Ruby] Play with the classic NGS simulation software wgsim ~Mysterious creature Pooh~

4 minute read

Introduction

Mysterious creature Poo It is a mysterious creature that may have existed since ancient times.

image.png Imagination of Pooh

What is wgsim?

wgsim is an NGS simulator created by giant Heng Li in the bioinfo world. It’s a 2011 software, but it works fine.

https://github.com/lh3/wgsim

Install

If you have Bioconda installed

conda install wgsim

~Mysterious creature Pooh~

image.png

The mysterious creature Pooh was born out of total whim. God creates creatures. (Evolution is outrageous.) Since you are a god, let’s create an organism with a sequence of 5000 bases, which is almost the same as a virus. The name of this creature is named “Poo”.

The reference array of Pooh is randomly generated. (I’m an endangered Ruby, but I can do the same with other programming languages.)

genome = Array.new(5000){ %w[A C T G].sample}
puts'>poo chr1'
genome.each_slice(80) do |l|
  puts l.join('')
end

When you run this script, you will get an array like this:

This is the mysterious creature **Pool’s reference genome! !! !! **

>poo chr1
AAAGCGCTGGATCATGCTTTAAGGCTTAGTTTTATGGCCCTAAACACGGAGACTGAAGCCGCTTGAGTCCGGTTTCAGAG
AGATTTGGTTTCACCCACAACTGCCAGCATATAATGGACCGTCGGCCGTACACGTTATGAGCCTGGCAGGGGTTACTACG
CGCCACGGGCTAGAGACCTAGTGGCGAAACGTGAAGGGCAATTATGGAATCCCGTACGCGCAATAAATGTCTACGCATTC
CCGGGGGGTGTTGTAGTCCCACGGCCCCGGACCTGCATTGGTCCAAATTCTCTTCGGTGCATATTAGCTTCATCGGAACC
GGCCTTTCTAATTCAGAACTGTTCGGTCATACACTATCTTCCAACACGCCGTGCACGAGTTGCACGGAGATCTAGTTGTT
ACTCGACGCGCGTACGACCGAGGCGGGGCTTAACACGCGAGTAGAAGCCCTGGGTTATCCACCACTTAATTTTGTAAAAA
GCGCTTTCAGTACTACACGGCTTGACCGTAGCTCTGTCGGACTTGTGACCCCTTCCGGCGAACATTAATTGTACCCCTTC
TTCTGACGGATTGGTTAGTAGTCCAAGCAGTCTTCACCCCAGGAGCACTTCTTCTTGCCGCGAGGCCATTCAGTCCAGTT
CGTGATGGCGGGGGTACCGATGCAAAAGCTCGTCTTTCTTGAGGCATCCACAGCCAGAGGCGGGGTATTTACGAGATTAA
GTTCTATCATTGCACTAGACAAAGCAAGTACAGCTCGGCAAAAAGAAACTCGTGGGCGCGTCTTGACCTGCAGCTAATAC
TGACTAGCATCCGTAGATATTACCGATCGTGCCGCACAAGGGAAATCTCGAGGCTGACAACGTCTTTGCCAAGGCAGTGG
CCGTGTCTTCTCGCCAGGGGTAAAGGTGCCGCTGGACTAGCGCTAACTTAGTGTCGTAATATCTCCATATGGACGATTTC
CCCGATCTTACCGGCTAGCCGCTGGTTGTATGCCGCGCCCAGGCCACCGAGTCTGACATGTTCCAAGGGCCTCATGTGAA
AGCGATGTGCTGATCCAGCATGGCCTTCAGCCTCACAAGACTGAGAGTTGGTGTAGAAACCCGGTTTGTTGGTAGTATGA
TGCTTAAAAATCGGTTTAATTTAGCTTGGTTAGATTTGAACGTGAAGAGCCTATGTCCAAAGACCAACTAACCTACTCAG
ACTAAGATTTACCGATTTCTCGTTTTGACTGAGATACGTAGTTAGCTATACGGGGCGTCGCCATCCTCTACCCCTTGAAT
ATGATTCGTAAGTGATGCCATTCTCGTGCTCGTGCGCTATGAATAGCAATGTCACCACAATAGTGCCCGGTGGCCGACAA
TCCGAGACCATGTCCGTCCTCATTCATGCGTAGTTAGTAGTGTATCGACTCACCTGGAGATGACCTCTATGTCCTGCTGC
GGTACGGCACGGCAGGATCTTTCTAGAAACATCGGCGATGCCCCAGAGGTTCTCTAAACTGTGTCATCTTTTTCTGCGGT
GCGTGGCTGGTCAAACCATTCAAATATATGAGAGTCATGCGCTAGCTCCCGGTAAAACCCCACAATAACGGGATTTGGGG
TTGAAAATGCGAGGCAGCGCACATGGATCAGGCCATTTAGGGTGTGATTGGTATGTTCCGAAAAAGCTGATGAGAACTAT
GCCCCGTTTACTTTGCCAGTAGAAAGAGTAGGAAGTAGTGTAACGCGGGCACCGGAAAACGACCGACACGATACTTTCCC
TTGTGCTGGAATTGTTTGATATATTAGGACGAGAATAGTTCTCTCCATTCCGCCGGCGTATAGACCAGGCCATGCCGTTC
TATTTCCAAGTTCCTTTGAAGGATCGGTGTTATTTCGGATGTCGTCGAGTAGTCGTGGGGCGTTAGGGTCTAAGCGACCA
GTTGTAGCTAATCGGGTGTCCTCCACTAACTGCTATTAGCATAGTACACCATGCTGCTTCAACTCATTTGGGATTTGCAC
TACTCTAGGCCAGAGATACAGATGCAAGCCGACTTAGCTGCACTAGCCAAGATTGAGCCCTAGGGCGTAGCGACCCTGAG
AATCCTGGCTAGCCCCCCTAGCCCCGGAAACTGCCCGTGTTGCTCTCGAAGCAGTTAGTAATCATCGTAATTTTGTCAAG
ATGGCCGTAGAGTGCAACCTGCAAGGCTGCAGGGTGTGTTCGGCATATTAATGCGGACAGTGCAGTGGACAAAGAGTCTC
GGAAGCGACGTGTATACAGTAGTGAGTAAATTGTTAAGCTTACTATCGCAGGAACTTCTGATCTGCCCGACCGTCGAGAT
CGAGTCTAGTAACAAGATTGCAGTCCTCAAGACTAGAAGTACCTAACTCCCTCAAATGCGAGTCCCAACAACTCATCCTC
ACCTCTGGAATCTCACCTCTAAGCATGACAGGCTAGGACCATTTAACATAGTTACAGTGACCTGATACGTAACGCCTCGA
AAGGAAAGGCGCGCCTAGGCAGCAATTGCAACGAACTCGGACGTCAGTACTCAACCAACCGAGTTATTTCAGCGTCGCGA
ATCTAATGCATCTTCCTGCTTAAGACGTCGTACTCCGGTGTCACCCTTCGACTAACTTAAAATCGCGCATTGGCAAGCAT
ACCTCTAATACCGGATCGCTCCCCTCAGGTTTACTTGATGAGCGCTGAACCCAACCATAATTTCGGCAAGGTAATCGCTT
ATCGTTCCACCCAGGCGTGAATAATTCATGCCCAATGCAGGTGACCGTATCTCAGGGAGATATAAGTGAAAGTTATGTGG
AAGGTATTGTAGTACATTTCAGGGTTCATCTAAGCTTCAGGACCACAGATGCAGAAGACGGGTTAAAGGACATGTCGGCA
CGAAGACATATCTTGTGGATCCCTATTCAGTGCCCGTGACTCAACAGATAGACTCATTTTCTCCTTTGGCGTGTCCGCAG
GTAAGATAACTTCAGTGTCCGAGCGGCGGAGTGGGTGTTCACGTCGTCGAGTGTCCATTTCGATTGGGGATGTTGGCTCG
CATACCACACGACCCTTTAGTGGACGGTCTAGTTTGTATACTCCCCGATTGGACTTTACCCGCGAAGCTAGCACTGTAAA
TATTGCCCCGGAGTGGCGGGATACGTTGGATATGCCTATACTGGACCGCAGCCTATTACATCAGATCGGCTTACCGCCAA
TCTTAGTGCTCAGCTACCCAGAGCTGTATCCATGGACATAGGCGCTATGACAAGCAGTACTGTAAAGCAGGGACAGCTTT
ACTATGTCGACAGCCGGCGTACGTCGTTGGGCTACCGCCGGCCTAGTTCCTAGTCAGTTCCGAGGCATAATCTACAAGTC
TAAGTAGGGCGCCATGATGTCACGAGACTTATAGGAGTCCACCGGCAAGTGCGTTGACGTAAATACATATGTGTATGCAG
CACTTTTGACAGGCCCGTAGATTTTGCACAAACAGCCGGGTGTGGTTCAAAAAGGGTGTTGCGGTCAGCCGCAGCAGCGC
AATTGTAACCAATCTACCTCCCGTTTGCCGACCGCTCGCGCAGATTCATGTAGGAACCGCTAGTTCCCCGTCCGCTACAT
TTCAATGGGCTTTTGCATGTCCTATCGCGATACCATCTCTGGTGGCGATGATATTTTATCAGGTTGCATTGGCCCAGTCC
ATCTTCGCGGCCGATTGGGCACGCCATTAGTTGCTCTCTGCGCAATCGCAGCTATACTAATGGTATGGTTGGGGAAATCA
TTGGAGTGCATTATTCCGCCCTCTTAACCTGTTGTGAACGGAACCGCAGTACGGTAGATTGCAAAGTCCAAGGTTTCTGG
ATTGTCTCTGAAAATGAAAGTTCTGCTCTGTCAAGTGATAAACGGAACGTCCTAGGAGGCGCCTACCGCGTTTAGAGAAG
TAACCCACGATATCTCGAACAAAAAATCAAAAGCGGTGTAGGTACAAAACGTATCGCTTACATCCCGTGAAGCTTAGTAGGACGCACCCTTCATCTCTCAGGGGGTCGCGTACACGATGTCCGTCTTTTATCGTGCCACGTGTACAAGTAGGGGTGCAGA
AATAACTGTACTAATGCAACGTTTCAGTTCTTCTGATGAACATCTCACACGGTTGTCCCCTCTAATATTAGGATTATCGG
CCTAGACGTGTTATGCACGCCGACCCCTTTCCCTAAGTGCTTCCGCAATGAAGTCGCGTATTTGTAACACACATTGGGTA
TGAATATAACTGTCCCTCCTGGAATAAGGAGTGCGCTTATCCTTTCGGACGTCAGAGCCCAGGGAGGAGACGCGTCAGAT
GCCATAGCCTATCGGGATCCGTCGTGCTAGTAGAACCAAAGTAACGGTGCTGTTAAATCGTTCCCTCAGCAGTCCGGGGT
TGACAGGTACAGTGTTTCCAACTACTAGAGGATCTCGCTGCTATAGAACGATACCACTGCGGGGGACACCGCCACTTTTC
GAGCTGTGGGGCCAGAAGTGCTGTCTGGATAATTCGGTAGGCATCAGCGGACAGCGGAGCTCGTAAGTGCCATGCACTTG
TAGCCCATTACGATGAGATCGGTTCCGCCAACTCCGTCAGCACAGGTATCCACGAGGAGTCACCGTCCGGCATGACGCAG
GACGGGGGCGGTGGAAATAATCGAGGCGGATACTCTTCATTTGTAGGCGCCAGCTACAGCGGTCCGGCGGGTCAATACCG
CATCCGCAGACATAATTGACCTCACTCAGCTCGCACGAACCTCCCTACTCTTTTCTACTCGAGGGTGCCTGATCAGAGAA
CGGAAGACGCTATAGGCAAGGCGGGCACGCCTGGAATACAGTCCCTTTCCGAATATTCACACCCGCTACCTATGTCGTAG
CGTTGCCGGGCCCGTTATCCAGACCGAACCCGTTGTCTGTTATAACATATTGAGTAATGGGCAGACAAACCCGTGTAGAG
TCCTCTCGGTTCCACATATCGAAGAGGGCGATCCGTACGA

Long.

Save it as poo.fa.

Running wgsim ~ Pooh’s whole genome sequence ~

By the way, the mysterious creature Poo was discovered in toothpaste for sheep. By the way, do sheep squirt? Now I don’t know the details. However, Pooh is reported to have been found in toothpaste that young sheep secretly used, so let’s believe it. It is said that the scientific name is pull-pull-pulls because it looks slick under a microscope. (I put it on right now)

Now, the person who discovered Pooh lost his interest in Pooh because of the mountains and valleys in his life, and he left him alone for 50 years. However, due to various accidental events such as the circulation of heaven and earth, the rise of the temperature on the equator, the chorus of the chairman of the town, the faint spring dad, the reconciliation of the chairman and the president, Pooh was sequenced with the whole genome. The time has come. Now, since you are a god, find the result of the genome sequence with wgsim.

Before running wgsim, let’s look at the options.

wgsim --help
Options: -e FLOAT base error rate [0.020]
         -d INT outer distance between the two ends [500]
         -s INT standard deviation [50]
         -N INT number of read pairs [1000000]
         -1 INT length of the first read [70]
         -2 INT length of the second read [70]
         -r FLOAT rate of mutations [0.0010]
         -R FLOAT fraction of indels [0.15]
         -X FLOAT probability an indel is extended [0.30]
         -S INT seed for random generator [0, use the current time]
         -A FLOAT discard if the fraction of ambiguous bases higher than FLOAT [0.05]
         -h haplotype mode

There are some options I don’t understand at all (such as what standard deviation is -s, or how -A or -h works), but I’ll study that later. (do not do)

The number of leads is 1000 here.

wgsim -N 1000 poo.fa poo_1.fq poo_2.fq

When executed, the mutation will be displayed. I think poo_1.fq and poo_2.fq have been generated. It’s easy.

poo 977 A W +
poo 1210 T W +
poo 1655 G K +
poo 1874 T K +
poo 2043 C S +
poo 2378 G T-

Let’s take a look at the contents of the generated fasta file.

head -n 12 poo_1.fq
@poo_418_960_2:0:0_1:0:0_0/1
TCGTCGAGTAGTCGTGGGGCGTTAGGGTCTAAGCGACCAGTTGTAGCTAATCGGGTGTCCTCCACTAACT
+
2222222222222222222222222222222222222222222222222222222222222222222222
@poo_624_1182_3:0:0_1:0:0_1/1
CTTAGGACATAGGCTCTTCACGTTCAAATCTAACCAAGCTAAATTAAACCGATTTTTAAGCATCATACTA
+
2222222222222222222222222222222222222222222222222222222222222222222222
@poo_2608_3166_3:0:0_0:0:0_2/1
GTCCAGTATAGGCATATCCAACGTATCCCGCCACTCCGGGGCAATATTTACAGTGCTAGCTTCGCGGGTA
+
2222222222222222222222222222222222222222222222222222222222222222222222

Well, I think the quality score of the second half of the normal lead will fall, but that is not the case. Let’s check the sequence quality with Fastqc.

image.png

You can see that the quality is constant (and low).

Generally, the quality should decrease as you move toward the back of the lead as shown in the figure below. (The image is from the official website of Fastqc https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)

Well, the problem there is old software, so it can’t be helped.

Mapping using BWA

Now, since it is a whole genome sequence, I will map it with bwa. What what. It’s a minimap? This is a world view from 10 years ago, and it’s an article played with wgsim, there is no long lead! !! I was upset. God must be careful not to be distraught.

First of all, create an index.

bwa index poo.fa

Pooh’s genome size is small, so I think execution will be completed in an instant. Let’s play bwa. (Bwa is also a tool made by Heng Li. It can be called God.)

bwa mem poo.fa poo_1.fq poo_2.fq >poo_wgs.sam

Well, let’s add SAM → BAM, BAM → sorted BAM, even indexes.

samtools view -bS poo_wgs.sam >poo_wgs.bam
samtools sort poo_wgs.bam poo_wgs.sorted
samtools index poo_wgs.sorted.bam

Let’s see the lead situation with IGV

In IGV, set poo.fa as the reference genome and open poo_was.sorted.bam.

Uooo! I’m done! I’m done! (Noisy)

igv_snapshot.png

Now, IGV actually has a convenient split view function. Let’s use this to observe the location of mutations.

image.png

Uooo! Older brother

The mutation has been properly detected.

in conclusion

By the way, I wrote an article about the mysterious creature ** Pooh ** in a whim, but it turned out that it was fun to write unexpectedly, and it was also a study. Pooh series, I may write if I feel like it again.