תופס כרומוזומים עבור hg38:
$ wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/*.fa.gz $ עבור fn ב- `ls * .fa.gz`; לעשות gunzip $ fn; נעשה
באמצעות מונה kmer ו- Python, כך תוכלו לחפש קמרים באורך 7 מכרומוזום chrY קוד>:
#! / usr / bin / env python import sysimport subprocessimport itertoolsk = 7chr = 'chrY'fastaFile ='% s.fa '% (chr) kmerCmd =' counter-counter --fasta --no-rc --k =% d% s '% (k, fastaFile) נסה: output = subprocess.check_output (kmerCmd, shell = True) result = {} עבור שורה ב- output.splitlines (): (כותרת, ספירה) = line.strip (). פיצול ('\ t') כותרת = כותרת [1:] kmers = dict ((key, int (val)) עבור (key, val) ב- [d.split ( ':') fo rd ב- counts.split ('')]) תוצאה [כותרת] = kmersexcept תת-תהליך. CalledProcessError כשגיאה: sys.stderr.write ("% s \ n"% (str (error))) kmers = תוצאה [chr] comp = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'} עבור kmerList ב- itertools.product ('ACGT', חזור על = k): kmerKey = '' .join (kmerList) kmerCompKey = '' .join (הפוך ([comp.get (b, b) עבור b ב- kmerList])) אם kmerKey לא ב- kmers ו- kmerCompKey לא ב- kmers: kmers [kmerKey] = 0 למפתח, val ב- מיון (kmers.iteritems (), key = lambda (key, val) :( val, key)): sys.stdout.write ("% s \ t% s \ n"% (key, val ))
סקריפט זה ידפיס קובץ טקסט המופרד באמצעות שתי עמודות לפלט רגיל, כאשר העמודה הראשונה היא ה- 7mer והעמודה השנייה היא ספירת ה- 7mer ההפוכה. משלימים מעל chrY (כולל ספירת אפס):
CGACGCG 20CGTCGCG 20CGCGATA 23TACGCGC 25 ... AAATAAA 33521TTTCTTT 34014GAATGGA 35361AATGGAA 36906TTTTTTTT 103093 em> וכרומוזומים לגנום הייחוס שלך.