[RUBY] Jouez avec le logiciel de simulation classique NGS wgsim ~ Mysterious creature Pooh ~
introduction
Créature mystérieuse ** Pooh ** C'est une ** créature mystérieuse qui peut avoir existé depuis les temps anciens **.
L'imagination de Poo
Qu'est-ce que wgsim?
wgsim est un simulateur NGS créé par Heng Li, un géant du monde de la bioinfo. C'est un logiciel de 2011, mais ça marche bien.
https://github.com/lh3/wgsim
Installation
Si vous avez installé Bioconda
conda install wgsim
~ Créature mystérieuse Pooh ~
La mystérieuse créature Pooh est née de caprices absolus. Dieu crée des créatures. (L'évolution est scandaleuse.) Puisque vous êtes un dieu, créons un organisme avec 5000 bases, qui est aussi séquencé qu'un petit virus. J'ai appelé cette créature "Poo".
La séquence de référence de Pooh est générée aléatoirement. (Je suis une ** secte Ruby ** en danger, mais vous pouvez faire de même avec d'autres langages de programmation.)
genome = Array.new(5000){ %w[A C T G].sample }
puts '>poo chr1'
genome.each_slice(80) do |l|
puts l.join('')
end
Si vous exécutez ce script, le tableau sera affiché comme ceci.
** Ceci est le génome de référence de la mystérieuse créature ** 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
Longue.
Enregistrez sous poo.fa.
Exécutez wgsim ~ Séquence complète du génome de Pooh ~
À propos, la mystérieuse créature ** Pooh ** a été découverte dans le dentifrice pour moutons. Au fait, les moutons tombent-ils malades? Eh bien, je ne connais pas les détails. Cependant, on rapporte que l'ourson a été trouvé dans le dentifrice secrètement utilisé par les jeunes moutons, alors croyons-le. On dit que le nom scientifique est pull-pull-pull car il semble dodu vu au microscope. (Je l'ai mis correctement maintenant)
Eh bien, la personne qui a découvert Pooh a perdu plus tard tout intérêt pour Pooh avec des montagnes et des vallées dans sa vie, et il est resté seul pendant les 50 années suivantes. Cependant, en raison de divers événements accidentels tels que le ciel et la terre, la montée de la température de l'équateur, le tremblement du président dans la ville, le léger conte printanier et la réconciliation entre le président et le président, Pooh a été séquencé dans tout le génome. Il est temps de venir. Maintenant, puisque vous êtes un dieu, utilisez wgsim pour trouver le résultat de la séquence du génome.
Avant d'exécuter wgsim, examinons les 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
Il y a certaines options dont je ne suis pas sûr (quel est l'écart type -s, comment fonctionne -A ou -h), mais j'étudierai cela plus tard. (ne pas faire)
Ici, le nombre de leads est fixé à 1000.
wgsim -N 1000 poo.fa poo_1.fq poo_2.fq
Une fois exécutée, la mutation sera affichée. Je pense que «poo_1.fq» et «poo_2.fq» ont été générés. C'est facile.
poo 977 A W +
poo 1210 T W +
poo 1655 G K +
poo 1874 T K +
poo 2043 C S +
poo 2378 G T -
Jetons un coup d'œil au contenu du fichier fasta généré.
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
Eh bien, je pense que le score de qualité baisse généralement dans la seconde moitié de l'avance, mais ce n'est pas le cas. Vérifions la qualité de la séquence avec Fastqc.
Vous pouvez voir que la qualité est constante (et faible).
En général, la qualité devrait diminuer à mesure que vous vous déplacez derrière la tête comme indiqué ci-dessous. (L'image provient du site officiel de Fastqc https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
Eh bien, il n'y a aucune aide, car il s'agit d'un ancien logiciel.
Cartographie à l'aide de BWA
Maintenant, puisqu'il s'agit d'une séquence complète du génome, nous allons la cartographier avec bwa. Quoi quoi. C'est la vision du monde il y a 10 ans, c'est un article à jouer avec wgsim, il n'y a pas de longue avance! !! Oups, je suis énervé. Dieu doit faire attention à ne pas être contrarié.
Tout d'abord, créez un index.
bwa index poo.fa
La taille du génome de Pooh est petite, donc je pense qu'il sera terminé en un instant. Appliquons bwa. (Bwa est également un outil fabriqué par Heng Li. Il peut être appelé un dieu.)
bwa mem poo.fa poo_1.fq poo_2.fq > poo_wgs.sam
Maintenant, ajoutons SAM → BAM, BAM → BAM trié, et même index.
samtools view -bS poo_wgs.sam > poo_wgs.bam
samtools sort poo_wgs.bam poo_wgs.sorted
samtools index poo_wgs.sorted.bam
Voyons comment Reed regarde IGV
Dans IGV, définissez poo.fa comme génome de référence et ouvrez poo_was.sorted.bam.
Uooo! C'est fait! J'ai fini, frère! (Bruyant)
Eh bien, IGV a en fait une fonction de vue fractionnée pratique. Utilisons ceci pour observer l'emplacement de la mutation.
Uooo! Frère (omis)
La mutation a été correctement détectée.
en conclusion
Eh bien, j'ai écrit un article sur la mystérieuse créature ** Pooh ** sur un coup de tête, mais il s'est avéré que c'était amusant à écrire et c'était aussi utile pour étudier. Série Pooh, peut-être que je l'écrirai à nouveau si j'en ai envie.