[RUBY] Spielen Sie mit der klassischen NGS-Simulationssoftware wgsim ~ Mysterious Kreatur Pooh ~

Einführung

Mysteriöse Kreatur ** Pooh ** Es ist eine ** mysteriöse Kreatur, die möglicherweise seit der Antike existiert hat **.

Poos Fantasie

Was ist wgsim?

wgsim ist ein NGS-Simulator, der von Heng Li, einem Riesen in der Bioinfo-Welt, entwickelt wurde. Es ist 2011 Software, aber es funktioniert gut.

https://github.com/lh3/wgsim

Installation

Wenn Sie Bioconda installiert haben

conda install wgsim

~ Geheimnisvolle Kreatur Pooh ~

Die mysteriöse Kreatur Pooh wurde aus Launen geboren. Gott erschafft Kreaturen. (Die Evolution ist empörend.) Da Sie ein Gott sind, erstellen wir einen Organismus mit 5000 Basen, der so sequenziert ist wie ein kleines Virus. Ich habe diese Kreatur "Poo" genannt.

Die Referenzsequenz von Pooh wird zufällig generiert. (Ich bin eine gefährdete ** Ruby-Sekte **, aber Sie können dasselbe mit anderen Programmiersprachen tun.)

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

Wenn Sie dieses Skript ausführen, wird das Array wie folgt ausgegeben.

** Dies ist das Referenzgenom der mysteriösen Kreatur ** Pooh **! !! !! ** **.

>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
TAACCCACGATATCTCGAACAAAAAATCAAAAGCGGTGTAGGTACAAAACGTATCGCTTACATCCCGTGAAGCTTAGTAG
GACGCACCCTTCATCTCTCAGGGGGTCGCGTACACGATGTCCGTCTTTTATCGTGCCACGTGTACAAGTAGGGGTGCAGA
AATAACTGTACTAATGCAACGTTTCAGTTCTTCTGATGAACATCTCACACGGTTGTCCCCTCTAATATTAGGATTATCGG
CCTAGACGTGTTATGCACGCCGACCCCTTTCCCTAAGTGCTTCCGCAATGAAGTCGCGTATTTGTAACACACATTGGGTA
TGAATATAACTGTCCCTCCTGGAATAAGGAGTGCGCTTATCCTTTCGGACGTCAGAGCCCAGGGAGGAGACGCGTCAGAT
GCCATAGCCTATCGGGATCCGTCGTGCTAGTAGAACCAAAGTAACGGTGCTGTTAAATCGTTCCCTCAGCAGTCCGGGGT
TGACAGGTACAGTGTTTCCAACTACTAGAGGATCTCGCTGCTATAGAACGATACCACTGCGGGGGACACCGCCACTTTTC
GAGCTGTGGGGCCAGAAGTGCTGTCTGGATAATTCGGTAGGCATCAGCGGACAGCGGAGCTCGTAAGTGCCATGCACTTG
TAGCCCATTACGATGAGATCGGTTCCGCCAACTCCGTCAGCACAGGTATCCACGAGGAGTCACCGTCCGGCATGACGCAG
GACGGGGGCGGTGGAAATAATCGAGGCGGATACTCTTCATTTGTAGGCGCCAGCTACAGCGGTCCGGCGGGTCAATACCG
CATCCGCAGACATAATTGACCTCACTCAGCTCGCACGAACCTCCCTACTCTTTTCTACTCGAGGGTGCCTGATCAGAGAA
CGGAAGACGCTATAGGCAAGGCGGGCACGCCTGGAATACAGTCCCTTTCCGAATATTCACACCCGCTACCTATGTCGTAG
CGTTGCCGGGCCCGTTATCCAGACCGAACCCGTTGTCTGTTATAACATATTGAGTAATGGGCAGACAAACCCGTGTAGAG
TCCTCTCGGTTCCACATATCGAAGAGGGCGATCCGTACGA

Lange.

Speichern Sie als poo.fa.

Führen Sie wgsim ~ Poohs gesamte Genomsequenz ~ aus

Übrigens wurde die mysteriöse Kreatur ** Pooh ** in Zahnpasta für Schafe entdeckt. Werden Schafe übrigens krank? Nun, ich kenne die Details nicht. Es wird jedoch berichtet, dass Pooh in der Zahnpasta gefunden wurde, die heimlich von jungen Schafen verwendet wird. Glauben wir es also. Der wissenschaftliche Name soll Pull-Pull-Pull sein, weil er unter dem Mikroskop prall aussieht. (Ich ziehe es jetzt richtig an)

Jetzt verlor die Person, die Pooh entdeckte, das Interesse an Pooh, der Berge und Täler in seinem Leben hatte, und verließ Pooh für die nächsten 50 Jahre. Aufgrund verschiedener zufälliger Ereignisse wie Himmel und Erde, des Temperaturanstiegs des Äquators, des Zitterns des Vorsitzenden in der Stadt, der schwachen Frühlingsgeschichte und der Versöhnung zwischen dem Vorsitzenden und dem Präsidenten wurde Pooh im gesamten Genom sequenziert. Es ist Zeit zu kommen. Nun, da Sie ein Gott sind, verwenden Sie wgsim, um das Ergebnis der Genomsequenz zu finden.

Bevor wir wgsim ausführen, werfen wir einen Blick auf die Optionen.

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

Es gibt einige Optionen, bei denen ich mir nicht sicher bin (welche Standardabweichung ist -s, wie -A oder -h funktioniert), aber ich werde das später untersuchen. (TU es nicht)

Hier ist die Anzahl der Ableitungen auf 1000 eingestellt.

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

Bei der Ausführung wird die Mutation angezeigt. Ich denke, dass poo_1.fq und poo_2.fq generiert wurden. Es ist einfach.

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

Werfen wir einen Blick auf den Inhalt der generierten Fasta-Datei.

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

Nun, ich denke, dass der Qualitätsfaktor normalerweise in der zweiten Hälfte der Führung abfällt, aber das ist nicht der Fall. Lassen Sie uns die Sequenzqualität mit Fastqc überprüfen.

Sie können sehen, dass die Qualität konstant (und niedrig) ist.

Im Allgemeinen sollte die Qualität abnehmen, wenn Sie sich wie unten gezeigt hinter die Leine bewegen. (Das Bild stammt von der offiziellen Website von Fastqc https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)

Nun, es gibt keine Hilfe dafür, weil es alte Software ist.

Mapping mit BWA

Da es sich nun um eine ganze Genomsequenz handelt, werden wir sie mit bwa abbilden. Was was. Was ist mit Minikarte? Dies ist die Weltanschauung vor 10 Jahren, es ist ein Artikel, um mit wgsim zu spielen, es gibt keinen langen Vorsprung! !! Ups, ich bin verärgert. Gott muss aufpassen, dass er nicht verärgert ist.

Erstellen Sie zunächst einen Index.

bwa index poo.fa

Die Genomgröße von Pooh ist klein, daher denke ich, dass sie in Kürze abgeschlossen sein wird. Wenden wir bwa an. (Bwa ist auch ein Werkzeug von Heng Li. Es kann als Gott bezeichnet werden.)

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

Fügen wir nun SAM → BAM, BAM → sortierte BAM und sogar Index hinzu.

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

Mal sehen, wie Reed IGV sieht

Setzen Sie in IGV "poo.fa" als Referenzgenom und öffnen Sie "poo_was.sorted.bam".

Uooo! Es ist fertig! Ich bin fertig, Bruder! (Laut)

Nun, IGV hat tatsächlich eine praktische Split-View-Funktion. Verwenden wir dies, um den Ort der Mutation zu beobachten.

Uooo! Bruder (weggelassen)

Die Mutation wurde korrekt erkannt.

abschließend

Nun, ich habe aus einer Laune heraus einen Artikel über die mysteriöse Kreatur ** Pooh ** geschrieben, aber es stellte sich heraus, dass das Schreiben Spaß machte und auch eine Lernerfahrung war. Pooh-Serie, vielleicht schreibe ich es noch einmal, wenn ich Lust dazu habe.