[RUBY] C'est un journal d'images des vacances d'été, où vous pouvez voir la séquence non spécifiée du génome humain.

introduction

Cet article est vraiment génial, alors ne le prenez pas vraiment.

Objectif

On dit que le programme du génome humain a révélé toute la séquence de l'ADN humain, mais il doit y avoir beaucoup de parties inconnues. Dessinez un graphique pour voir où vous ne savez pas.

Télécharger le génome humain fasta

Vous pouvez télécharger la séquence du génome humain depuis Gencode. https://www.gencodegenes.org/human/

Ici, téléchargez TOUS luxueusement.

wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_35/GRCh38.p13.genome.fa.gz

Vous pouvez le décompresser, mais vous pouvez le voir dans une certaine mesure en utilisant la commande zcat.

Tout d'abord, affichons head.

zcat GRCh38.p13.genome.fa.gz | head

>chr1 1
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

Oui, il est sorti immédiatement. C'est une montagne de N. En d'autres termes, la partie où la séquence ne peut pas être spécifiée vient du début. Ensuite, vérifions les chromosomes pour le moment.

zcat GRCh38.p13.genome.fa.gz | grep "^>"

Eh bien, il sort en grandes quantités comme ça.

Ensuite, vérifions si le tableau se compose vraiment de 5 types de chaînes qui ne sont que TCGAN. Je suis sûr qu'il existe un moyen d'utiliser les outils de ligne de commande, mais ici je vais écrire un petit programme en utilisant le langage Crystal. La raison du langage Crystal est qu'il est explosif.

count.cr

counter = Hash(Char, Int32).new(0)

while l = gets
    l.chomp.each_char do |c|
        counter[c] += 1
    end
end

p counter

Oui, c'est juste du code simple pour quiconque comprend Ruby. Le seul point est que vous avez spécifié le type lorsque vous avez généré Hash. Tout le reste est le même que Ruby. C'est Crystal. Construire.

crystal build counter.cr --release

Si vous incluez ce temps de construction, la vitesse d'exécution peut ne pas être très différente de celle de Julia etc., mais je ne m'inquiéterai pas à ce sujet. Comptons le nombre de TCGAN

zcat GRCh38.p13.genome.fa.gz | grep -v "^>" | ./count

La ligne qui définit le chromosome au début de la ligne est omise avec grep -v" ^> ". Le résultat est

{'N' => 161331703, 'T' => 916993511, 'A' => 914265135, 'C' => 635937481, 'G' => 638590158}

est devenu. Certes, il a été confirmé qu'il n'y a pas de chaînes autres que TCGA et N. Ensuite, comptons la distribution continue de N. En d'autres termes, découvrons combien il y a de N caractères consécutifs.

nseq.cr

temp = 0 

while l = gets
    l.chomp.each_char do |c|
        if c == 'N'
             temp += 1
        elsif temp != 0
            puts temp
            temp = 0
        end
    end
end
puts temp if temp != 0

Construire.

crystal build nseq.cr --release

Voyons d'abord combien il y a de N points consécutifs.

zcat GRCh38.p13.genome.fa.gz | grep -v "^>" | ./nseq | wc -l

1234

Vous pouvez voir qu'il n'y en a pas tellement. Si vous les affichez côte à côte dans l'ordre du plus long, cela ressemblera à ceci.

Je trouve étrange que 1 le soit, et j'ai l'impression qu'il y a beaucoup de centaines et de 50 000.

Sujet principal

Maintenant, écrivez un programme qui calcule le rapport de N, AT et CG pour 10 000 caractères.

n2

tcgan = Hash(Char, Int32).new(0)

target = ARGV[0]
chr = ""
loc = 1
flag = false

while l = gets
  if l.starts_with?(">")
    exit if flag
    if l == target
      puts "loc\tN\tAT\tCG"
      flag = true
    end
    next
  end
  if flag
    l.chomp.each_char do |c|
      tcgan[c] += 1
      loc += 1
      if loc % 10000 == 0
        total = tcgan.values.sum.to_f
        ta = (tcgan['A'] + tcgan['T']) / total
        cg = (tcgan['G'] + tcgan['C']) / total
        n  = tcgan['N'] / total
        puts "#{loc}\t#{n}\t#{ta}\t#{cg}"
        tcgan = Hash(Char, Int32).new(0)
      end
    end
  end
end

Exécutez-le comme un essai.

zcat GRCh38.p13.genome.fa.gz | ./n2 ">chr1 1" | head

loc	N	AT	CG
10000	1.0	0.0	0.0
20000	0.0001	0.4076	0.5923
30000	0.0	0.4826	0.5174
40000	0.0	0.5288	0.4712
50000	0.0	0.6439	0.3561
60000	0.0	0.6346	0.3654
70000	0.0	0.6669	0.3331
80000	0.0	0.6199	0.3801
90000	0.0	0.6294	0.3706

Il semble que cela fonctionne bien. La création d'un programme avec le même comportement dans Ruby prend beaucoup de temps à exécuter, mais Crystal est très rapide. Lancez ceci pour uplot sur la ligne de commande.

uplot est un outil Ruby que je crée personnellement et qui vous permet d'afficher des graphiques sur votre terminal en utilisant UnicodePlots.rb.

À partir de là, je fais quelque chose d'assez déroutant pour dessiner un graphique sur le terminal, prendre une capture d'écran et la coller dans Qiita.

chr1

chr2

chr3

chr4

chr5

chr6

chr7

chr8

chr9

chr10

chr11

chr12

chr13

chr14

chr 15

chr 16

chr 17

chr 18

chr 19

chr 20

chr 21

Ajout de dispersion

chr 22

Ajout de dispersion

chr X

chr Y

Selon Wikipedia Il semble que la région de l'hétérochromatine est principalement non déchiffrée et que cette région est inactive. Je pense que cela signifie qu'il est rarement transféré. En d'autres termes, cela peut ne pas être biologiquement significatif.

Cependant, quand j'essaye ceci, il semble qu'il y ait beaucoup de parties où l'arrangement est inconnu.

Si vous le connaissez, n'hésitez pas à y participer.

C'est tout pour cet article.