U kunt de Jellyfish -software gebruiken om de k-mer-profielen tot een lengte van 31 te berekenen.

Op basis van de instructies in de gebruiker guide:

Het basiscommando om alle k-mers te tellen is als volgt:

  kwallen tellen -m 21 -s 100M -t 10 -C reads.fasta  

Om het histogram van de k-mer occurrences te berekenen, gebruik je het histo-subcommando (zie paragraaf 3.1):

  jellyfish histo mer_counts.jf 

Om de tellingen van een bepaalde k-mer op te vragen, gebruik je de query-subopdracht (zie sectie 3.3):

  kwallenquery mer_counts.jf AACGTTG  

Om alle tellingen voor alle k-mers in het bestand uit te voeren, gebruik je het dump subcommando (zie paragraaf 3.2):

  jellyfish dump mer_counts.jf > mer_counts_dumps.fa  

Je kunt ook R. gebruiken. Ik geef je een voorbeeld van alleen chr1 en alleen kmer = 4.

  bibliotheek (BSgenome.Hsapiens.UCSC.hg38) bibliotheek (Biostrings) genoom <- BSgenome.Hsapiens.UCSC.hg38kmers <- oligonucleotideFrequency (genoom $ chr1, 4) kmersm <- as.matrix (kmers) m [order (m),]  

Chromosomen grijpen voor hg38:

  $ wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/*.fa.gz $ voor fn in `ls * .fa.gz`; doe gunzip $ fn; done  

Via kmer-counter en Python, hier is hoe te zoeken naar kmers met lengte 7 vanaf chromosoom chrY :

  #! / usr / bin / env python import sysimport subprocessimport itertoolsk = 7chr = 'chrY'fastaFile ='% s.fa '% (chr) kmerCmd =' kmer-counter --fasta --no-rc --k =% d% s '% (k, fastaFile) probeer: output = subprocess.check_output (kmerCmd, shell = True) resultaat = {} voor regel in output.splitlines (): (header, counts) = line.strip (). split ('\ t') header = header [1:] kmers = dict ((key, int (val)) voor (key, val) in [d.split ( ':') voor rd in counts.split ('')]) resultaat [header] = kmers behalve subprocess.CalledProcessError als fout: sys.stderr.write ("% s \ n"% (str (error))) kmers = resultaat [chr] comp = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'} voor kmerList in itertools.product ('ACGT', repeat = k): kmerKey = '' .join (kmerList) kmerCompKey = '' .join (omgekeerd ([comp.get (b, b) voor b in kmerList])) als kmerKey niet in kmers en kmerCompKey niet in kmers: kmers [kmerKey] = 0voor sleutel, waarde gesorteerd (kmers.iteritems (), key = lambda (key, val) :( val, key)): sys.stdout.write ("% s \ t% s \ n"% (key, val ))  

Dit script zal een tweekoloms door tabs gescheiden tekstbestand afdrukken naar standaarduitvoer, waarbij de eerste kolom de 7mer is en de tweede kolom de telling van die 7mer en het omgekeerde aanvulling op chrY (inclusief nultellingen):

  CGACGCG 20CGTCGCG 20CGCGATA 23TACGCGC 25 ... AAATAAA 33521TTTCTTT 34014GAATGGA 35361AATGGAA 36906TTTTTTT 103093  

Met een paar aanpassingen kan dit script over een bereik van em> en chromosomen voor uw referentiegenoom.