Βιοπληροφορική είναι ένα διεπιστημονικό πεδίο που αναπτύσσει μεθόδους και εργαλεία λογισμικού για την ανάλυση και κατανόηση των βιολογικών δεδομένων.

Με πιο απλά λόγια, μπορείτε απλά να το σκεφτείτε ως επιστημονικά δεδομένα για τη βιολογία.

Μεταξύ των πολλών τύπων βιολογικών δεδομένων, τα δεδομένα γονιδιωματικής είναι ένα από τα πιο ευρέως αναλυμένα. Ειδικά με την ταχεία πρόοδο των τεχνολογιών αλληλουχίας DNA επόμενης γενιάς (NGS), ο όγκος των δεδομένων γονιδιωματικής αυξάνεται εκθετικά. Σύμφωνα με Stephens, Zachary D et αϊ. , η απόκτηση δεδομένων γονιδιωματικής είναι σε κλίμακα exabyte ανά έτος.

Η SciPy συγκεντρώνει πολλές μονάδες Python για επιστημονικούς υπολογιστές, οι οποίες είναι ιδανικές για πολλές βιοπληροφορικές ανάγκες. Τιτίβισμα

Σε αυτήν την ανάρτηση, επιδεικνύω ένα παράδειγμα ανάλυσης ενός αρχείου GFF3 για το ανθρώπινο γονιδίωμα με το SciPy Stack. Generic Feature Format Έκδοση 3 (GFF3) είναι η τρέχουσα τυπική μορφή αρχείου κειμένου για την αποθήκευση γονιδιωματικών χαρακτηριστικών. Συγκεκριμένα, σε αυτήν την ανάρτηση θα μάθετε πώς να χρησιμοποιείτε τη στοίβα SciPy για να απαντήσετε στις ακόλουθες ερωτήσεις σχετικά με το ανθρώπινο γονιδίωμα:

1. Πόσο από το γονιδίωμα είναι ατελές;
2. Πόσα γονίδια υπάρχουν στο γονιδίωμα;
3. Πόσο καιρό είναι ένα τυπικό γονίδιο;
4. Πώς μοιάζει η κατανομή των γονιδίων μεταξύ των χρωμοσωμάτων;

Μπορείτε να κατεβάσετε το πιο πρόσφατο αρχείο GFF3 για το ανθρώπινο γονιδίωμα εδώ . ο ΔΙΑΒΑΣΤΕ Το αρχείο που περιλαμβάνεται σε αυτόν τον κατάλογο παρέχει μια σύντομη περιγραφή αυτής της μορφής δεδομένων και βρίσκεται πιο λεπτομερής περιγραφή εδώ .

Θα το χρησιμοποιησουμε Πάντες , ένα σημαντικό στοιχείο της στοίβας SciPy που παρέχει γρήγορες, ευέλικτες και εκφραστικές δομές δεδομένων, για χειρισμό και κατανόηση του αρχείου GFF3.

Ρύθμιση

Πρώτα τα πράγματα πρώτα, ας ρυθμίσουμε ένα εικονικό περιβάλλον με εγκατεστημένη τη στοίβα SciPy. Αυτή η διαδικασία μπορεί να είναι χρονοβόρα εάν δημιουργηθεί από την πηγή με μη αυτόματο τρόπο, καθώς η στοίβα περιλαμβάνει πολλά πακέτα - μερικά από τα οποία εξαρτώνται από εξωτερικό κώδικα FORTRAN ή C. Εδώ, προτείνω τη χρήση Μινικόντα , γεγονός που καθιστά την εγκατάσταση πολύ εύκολη.

wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh bash Miniconda3-latest-Linux-x86_64.sh -b

Το -b flag στη γραμμή bash του λέει να εκτελεστεί σε batch mode. Αφού χρησιμοποιηθούν οι παραπάνω εντολές για την επιτυχή εγκατάσταση του Miniconda, ξεκινήστε ένα νέο εικονικό περιβάλλον για τη γονιδιωματική και, στη συνέχεια, εγκαταστήστε τη στοίβα SciPy.

mkdir -p genomics cd genomics conda create -p venv ipython matplotlib pandas

Σημειώστε ότι έχουμε καθορίσει μόνο τα 3 πακέτα που θα χρησιμοποιήσουμε σε αυτήν την ανάρτηση. Εάν θέλετε όλα τα πακέτα που αναφέρονται στη στοίβα SciPy, απλώς προσθέστε τα στο τέλος του conda create εντολή.

Εάν δεν είστε σίγουροι για το ακριβές όνομα ενός πακέτου, δοκιμάστε conda search. Ας ενεργοποιήσουμε το εικονικό περιβάλλον και ξεκινήστε το IPython.

source activate venv/ ipython

IPython είναι μια πολύ πιο ισχυρή αντικατάσταση της προεπιλογής Διεπαφή διερμηνέα Python , οπότε ό, τι κάνατε στον προεπιλεγμένο διερμηνέα python μπορεί επίσης να γίνει στο IPython. Συνιστώ ανεπιφύλακτα σε κάθε προγραμματιστή Python, που δεν έχει χρησιμοποιήσει ακόμα το IPython, να το δοκιμάσει.

Κατεβάστε το αρχείο σχολιασμού

Με την ολοκλήρωση της ρύθμισής μας, ας κατεβάσουμε το αρχείο σχολιασμού ανθρώπινου γονιδιώματος σε μορφή GFF3.

Είναι περίπου 37 MB, ένα πολύ μικρό αρχείο σε σύγκριση με το περιεχόμενο πληροφοριών ενός ανθρώπινου γονιδιώματος, το οποίο είναι περίπου 3 GB σε απλό κείμενο. Αυτό συμβαίνει επειδή το αρχείο GFF3 περιέχει μόνο τον σχολιασμό των ακολουθιών, ενώ τα δεδομένα ακολουθίας συνήθως αποθηκεύονται σε άλλη μορφή αρχείου που ονομάζεται FASTA . Αν σας ενδιαφέρει, μπορείτε να κατεβάσετε το FASTA εδώ , αλλά δεν θα χρησιμοποιήσουμε τα δεδομένα ακολουθίας σε αυτό το σεμινάριο.

!wget ftp://ftp.ensembl.org/pub/release-85/gff3/homo_sapiens/Homo_sapiens.GRCh38.85.gff3.gz

Το πρόθεμα ! λέει στο IPython ότι αυτή είναι μια εντολή shell αντί για μια εντολή Python. Ωστόσο, το IPython μπορεί επίσης να επεξεργαστεί μερικές συχνά χρησιμοποιούμενες εντολές κελύφους όπως ls, pwd, rm, mkdir, rmdir ακόμη και χωρίς πρόθεμα !.

Ρίξτε μια ματιά στο κεφάλι του αρχείου GFF, θα δείτε πολλές γραμμές μεταδεδομένων / πραγμάτων / οδηγιών που ξεκινούν με ## ή #!.

Σύμφωνα με την ΔΙΑΒΑΣΤΕ , ## σημαίνει ότι τα μεταδεδομένα είναι σταθερά, ενώ #! σημαίνει ότι είναι πειραματικό.

Αργότερα θα δείτε επίσης ###, η οποία είναι μια άλλη οδηγία με ακόμη πιο λεπτή έννοια που βασίζεται στο προσδιορισμός .

Τα ανθρώπινα αναγνώσιμα σχόλια υποτίθεται ότι ακολουθούν ένα μόνο #. Για απλότητα, θα αντιμετωπίσουμε όλες τις γραμμές ξεκινώντας από # ως σχόλια και απλώς αγνοήστε τα κατά την ανάλυσή μας.

##gff-version 3 ##sequence-region 1 1 248956422 ##sequence-region 10 1 133797422 ##sequence-region 11 1 135086622 ##sequence-region 12 1 133275309 ... ##sequence-region MT 1 16569 ##sequence-region X 1 156040895 ##sequence-region Y 2781480 56887902 #!genome-build GRCh38.p7 #!genome-version GRCh38 #!genome-date 2013-12 #!genome-build-accession NCBI:GCA_000001405.22 #!genebuild-last-updated 2016-06

Η πρώτη γραμμή δείχνει ότι η έκδοση της μορφής GFF που χρησιμοποιείται σε αυτό το αρχείο είναι 3.

Ακολουθούν περιλήψεις όλων των περιοχών ακολουθίας. Όπως θα δούμε αργότερα, τέτοιες πληροφορίες μπορούν επίσης να βρεθούν στο κύριο μέρος του αρχείου.

Οι γραμμές που ξεκινούν με #! εμφανίζει πληροφορίες σχετικά με τη συγκεκριμένη έκδοση του γονιδιώματος, GRCh38.p7, στην οποία ισχύει αυτό το αρχείο σχολιασμού.

Κοινοπραξία αναφοράς γονιδιώματος (GCR) είναι μια διεθνής κοινοπραξία, η οποία επιβλέπει ενημερώσεις και βελτιώσεις σε πολλά συγκροτήματα γονιδιώματος αναφοράς, συμπεριλαμβανομένων εκείνων για ανθρώπους, ποντίκια, ζέβρα και κοτόπουλα.

Σάρωση μέσω αυτού του αρχείου, ακολουθούν οι πρώτες γραμμές σχολιασμού.

Παραδείγματα προϋπολογισμού κεφαλαίου με λύσεις

1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_000001.11 ### 1 . biological_region 10469 11240 1.3e+03 . . external_name=oe %3D 0.79;logic_name=cpg 1 . biological_region 10650 10657 0.999 + . logic_name=eponine 1 . biological_region 10655 10657 0.999 - . logic_name=eponine 1 . biological_region 10678 10687 0.999 + . logic_name=eponine 1 . biological_region 10681 10688 0.999 - . logic_name=eponine ...

Οι στήλες είναι seqid , πηγή , τύπος , αρχή , τέλος , σκορ , Παραλία , φάση , γνωρίσματα . Μερικά από αυτά είναι πολύ εύκολα κατανοητά. Πάρτε την πρώτη γραμμή ως παράδειγμα:

1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_000001.11

Αυτός είναι ο σχολιασμός του πρώτου χρωμοσώματος με seqid 1, που ξεκινά από την πρώτη βάση έως την 24.895.622η βάση.

Με άλλα λόγια, το πρώτο χρωμόσωμα έχει μήκος περίπου 25 εκατομμύρια βάσεις.

Η ανάλυσή μας δεν θα χρειαστεί πληροφορίες από τις τρεις στήλες με τιμή . (δηλαδή σκορ, σκέλος και φάση), έτσι μπορούμε απλώς να τους αγνοήσουμε προς το παρόν.

Η τελευταία στήλη χαρακτηριστικών λέει ότι το Chromosome 1 έχει επίσης τρία ψευδώνυμα ονόματα, δηλαδή CM000663.2, chr1 και NC_000001.11. Βασικά είναι αυτό που μοιάζει ένα αρχείο GFF3, αλλά δεν θα τα ελέγξουμε κάθε γραμμή, οπότε ήρθε η ώρα να φορτώσετε ολόκληρο το αρχείο στο Pandas.

Το Pandas είναι κατάλληλο για την αντιμετώπιση της μορφής GFF3 επειδή είναι ένα αρχείο οριοθετημένο με καρτέλες και το Pandas έχει πολύ καλή υποστήριξη για την ανάγνωση αρχείων τύπου CSV.

Σημειώστε ότι μια εξαίρεση στη μορφή οριοθετημένης καρτέλας είναι όταν το GFF3 περιέχει ##FASTA .

Σύμφωνα με την προσδιορισμός , ##FASTA υποδεικνύει το τέλος ενός τμήματος σχολιασμού, το οποίο θα ακολουθηθεί με μία ή περισσότερες ακολουθίες σε μορφή FASTA (χωρίς καρτέλα). Αλλά αυτό δεν ισχύει για το αρχείο GFF3 που πρόκειται να αναλύσουμε.

In [1]: import pandas as pd In [2]: pd.__version__ Out[2]: '0.18.1' In [3]: col_names = ['seqid', 'source', 'type', 'start', 'end', 'score', 'strand', 'phase', 'attributes'] Out[3]: df = pd.read_csv('Homo_sapiens.GRCh38.85.gff3.gz', compression='gzip', sep=' ', comment='#', low_memory=False, header=None, names=col_names)

Η τελευταία γραμμή παραπάνω φορτώνει ολόκληρο το αρχείο GFF3 με pandas.read_csv μέθοδος.

Δεδομένου ότι δεν είναι ένα τυπικό αρχείο CSV, πρέπει να προσαρμόσουμε λίγο την κλήση.

Αρχικά, ενημερώνουμε τους Pandas για τη μη διαθεσιμότητα πληροφοριών κεφαλίδας στο GFF3 με header=None και στη συνέχεια καθορίζουμε το ακριβές όνομα για κάθε στήλη με names=col_names.

Εάν το names Το όρισμα δεν έχει καθοριστεί, τα Pandas θα χρησιμοποιούν αυξητικούς αριθμούς ως ονόματα για κάθε στήλη.

sep=' ' λέει στους Pandas ότι οι στήλες διαχωρίζονται με καρτέλες αντί για διαχωρισμένες με κόμμα. Η τιμή σε sep= μπορεί στην πραγματικότητα να είναι μια κανονική έκφραση (regex). Αυτό γίνεται πρακτικό εάν το αρχείο χρησιμοποιεί διαφορετικά διαχωριστικά για κάθε στήλη (ω ναι, αυτό συμβαίνει). comment='#' σημαίνει γραμμές που ξεκινούν με # θεωρούνται σχόλια και θα αγνοηθούν.

compression='gzip' λέει στον Pandas ότι το αρχείο εισαγωγής είναι συμπιεσμένο με gzip.

Επιπλέον, pandas.read_csv έχει ένα πλούσιο σύνολο παραμέτρων που επιτρέπει την ανάγνωση διαφορετικών τύπων μορφών αρχείων τύπου CSV.

Ο τύπος της επιστρεφόμενης τιμής είναι a DataFrame, η οποία είναι η πιο σημαντική δομή δεδομένων στο Pandas, που χρησιμοποιείται για την αναπαράσταση 2D δεδομένων.

Το Pandas έχει επίσης ένα Series και Panel δομή δεδομένων για δεδομένα 1D και 3D, αντίστοιχα. Ανατρέξτε στο τεκμηρίωση για μια εισαγωγή στις δομές δεδομένων του Pandas.

Ας ρίξουμε μια ματιά στις πρώτες καταχωρήσεις με .head μέθοδος.

In [18]: df.head() Out[18]: seqid source type start end score strand phase attributes 0 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_00000... 1 1 . biological_region 10469 11240 1.3e+03 . . external_name=oe %3D 0.79;logic_name=cpg 2 1 . biological_region 10650 10657 0.999 + . logic_name=eponine 3 1 . biological_region 10655 10657 0.999 - . logic_name=eponine 4 1 . biological_region 10678 10687 0.999 + . logic_name=eponine

Η έξοδος διαμορφώνεται όμορφα σε μορφή πίνακα με μεγαλύτερες συμβολοσειρές στη στήλη χαρακτηριστικών αντικαθίστανται εν μέρει με ....

Μπορείτε να ρυθμίσετε τα Pandas να μην παραλείπουν μεγάλες χορδές με pd.set_option('display.max_colwidth', -1). Επιπλέον, το Pandas έχει πολλά επιλογές που μπορεί να προσαρμοστεί.

Στη συνέχεια, ας πάρουμε μερικές βασικές πληροφορίες σχετικά με αυτό το πλαίσιο δεδομένων με το .info μέθοδος.

In [20]: df.info() RangeIndex: 2601849 entries, 0 to 2601848 Data columns (total 9 columns): seqid object source object type object start int64 end int64 score object strand object phase object attributes object dtypes: int64(2), object(7) memory usage: 178.7+ MB

Αυτό δείχνει ότι ο GFF3 έχει 2.601.848 σχολιασμένες γραμμές και κάθε γραμμή έχει εννέα στήλες.

Για κάθε στήλη, δείχνει επίσης τους τύπους δεδομένων τους.

Αυτή η αρχή και το τέλος είναι τύπου int64, ακέραιοι που αντιπροσωπεύουν θέσεις στο γονιδίωμα.

Οι άλλες στήλες είναι όλες τύπου object, πράγμα που σημαίνει ότι οι τιμές τους αποτελούνται από ένα μείγμα ακέραιων αριθμών, πλωτών και συμβολοσειρών.

Το μέγεθος όλων των πληροφοριών είναι περίπου 178,7+ MB αποθηκευμένα στη μνήμη. Αυτό αποδεικνύεται πιο συμπαγές από το ασυμπίεστο αρχείο, το οποίο θα είναι περίπου 402 MB. Μια γρήγορη επαλήθευση εμφανίζεται παρακάτω.

gunzip -c Homo_sapiens.GRCh38.85.gff3.gz > /tmp/tmp.gff3 && du -s /tmp/tmp.gff3 402M /tmp/tmp.gff3

Από μια προβολή υψηλού επιπέδου, έχουμε φορτώσει ολόκληρο το αρχείο GFF3 σε ένα αντικείμενο DataFrame στην Python και όλες οι ακόλουθες αναλύσεις μας θα βασίζονται σε αυτό το μεμονωμένο αντικείμενο.

Τώρα, ας δούμε τι είναι η πρώτη στήλη seqid.

In [29]: df.seqid.unique() Out[29]: array(['1', '10', '11', '12', '13', '14', '15', '16', '17', '18', '19', '2', '20', '21', '22', '3', '4', '5', '6', '7', '8', '9', 'GL000008.2', 'GL000009.2', 'GL000194.1', 'GL000195.1', ... 'KI270757.1', 'MT', 'X', 'Y'], dtype=object) In [30]: df.seqid.unique().shape Out[30]: (194,)

df.seqid είναι ένας τρόπος πρόσβασης στα δεδομένα στηλών από ένα πλαίσιο δεδομένων. Ένας άλλος τρόπος είναι df['seqid'], η οποία είναι πιο γενική σύνταξη, επειδή εάν το όνομα της στήλης είναι δεσμευμένη λέξη-κλειδί Python (π.χ. class) ή περιέχει . ή χαρακτήρα διαστήματος, ο πρώτος τρόπος (df.seqid) δεν θα λειτουργήσει.

Η έξοδος δείχνει ότι υπάρχουν 194 μοναδικά seqids, τα οποία περιλαμβάνουν Chromosomes 1 έως 22, X, Y και DNA μιτοχονδρίων (MT), καθώς και 169 άλλες σειρές.

Οι σειρές που ξεκινούν με ΚΙ και GL είναι αλληλουχίες DNA - ή ικριώματα - στο γονιδίωμα που δεν έχουν συναρμολογηθεί επιτυχώς στα χρωμοσώματα.

Για όσους δεν είναι εξοικειωμένοι με τη γονιδιωματική, αυτό είναι σημαντικό.

Αν και το πρώτο σχέδιο ανθρώπινου γονιδιώματος βγήκε πριν από περισσότερα από 15 χρόνια, το τρέχον ανθρώπινο γονιδίωμα εξακολουθεί να είναι ατελές. Η δυσκολία στη συναρμολόγηση αυτών των ακολουθιών οφείλεται σε μεγάλο βαθμό πολύπλοκες επαναλαμβανόμενες περιοχές στο γονιδίωμα .

Στη συνέχεια, ας ρίξουμε μια ματιά στη στήλη προέλευσης.

ο ΔΙΑΒΑΣΤΕ αναφέρει ότι η πηγή είναι ένας προσδιοριστής ελεύθερου κειμένου που προορίζεται να περιγράψει τον αλγόριθμο ή τη διαδικασία λειτουργίας που δημιούργησε αυτήν τη δυνατότητα.

In [66]: df.source.value_counts() Out[66]: havana 1441093 ensembl_havana 745065 ensembl 228212 . 182510 mirbase 4701 GRCh38 194 insdc 74

Αυτό είναι ένα παράδειγμα της χρήσης του value_counts μέθοδος, η οποία είναι εξαιρετικά χρήσιμη για μια γρήγορη καταμέτρηση κατηγορηματικών μεταβλητών.

Από το αποτέλεσμα, μπορούμε να δούμε ότι υπάρχουν επτά πιθανές τιμές για αυτήν τη στήλη και η πλειονότητα των καταχωρήσεων στο αρχείο GFF3 προέρχεται από την Αβάνα, το ensembl και το ensembl_havana.

Μπορείτε να μάθετε περισσότερα σχετικά με το τι σημαίνουν αυτές οι πηγές και τις σχέσεις μεταξύ τους αυτήν την ανάρτηση .

Για να διατηρήσουμε τα πράγματα απλά, θα επικεντρωθούμε σε καταχωρήσεις από πηγές GRCh38, Αβάνα, ensembl και ensembl_havan.a.

Πόσο μεγάλο μέρος του γονιδιώματος είναι ελλιπές;

Οι πληροφορίες για κάθε ολόκληρο χρωμόσωμα βρίσκονται στις καταχωρίσεις από την πηγή GRCh38, οπότε ας πρώτα φιλτράρουμε τα υπόλοιπα και αντιστοιχίζουμε το φιλτραρισμένο αποτέλεσμα σε μια νέα μεταβλητή gdf.

In [70]: gdf = df[df.source == 'GRCh38'] In [87]: gdf.shape Out[87]: (194, 9) In [84]: gdf.sample(10) Out[84]: seqid source type start end score strand phase attributes 2511585 KI270708.1 GRCh38 supercontig 1 127682 . . . ID=supercontig:KI270708.1;Alias=chr1_KI270708v... 2510840 GL000208.1 GRCh38 supercontig 1 92689 . . . ID=supercontig:GL000208.1;Alias=chr5_GL000208v... 990810 17 GRCh38 chromosome 1 83257441 . . . ID=chromosome:17;Alias=CM000679.2,chr17,NC_000... 2511481 KI270373.1 GRCh38 supercontig 1 1451 . . . ID=supercontig:KI270373.1;Alias=chrUn_KI270373... 2511490 KI270384.1 GRCh38 supercontig 1 1658 . . . ID=supercontig:KI270384.1;Alias=chrUn_KI270384... 2080148 6 GRCh38 chromosome 1 170805979 . . . ID=chromosome:6;Alias=CM000668.2,chr6,NC_00000... 2511504 KI270412.1 GRCh38 supercontig 1 1179 . . . ID=supercontig:KI270412.1;Alias=chrUn_KI270412... 1201561 19 GRCh38 chromosome 1 58617616 . . . ID=chromosome:19;Alias=CM000681.2,chr19,NC_000... 2511474 KI270340.1 GRCh38 supercontig 1 1428 . . . ID=supercontig:KI270340.1;Alias=chrUn_KI270340... 2594560 Y GRCh38 chromosome 2781480 56887902 . . . ID=chromosome:Y;Alias=CM000686.2,chrY,NC_00002...

Το φιλτράρισμα είναι εύκολο στο Pandas.

Εάν ελέγξετε την τιμή που αξιολογήθηκε από την έκφραση df.source == 'GRCh38', είναι μια σειρά True και False τιμές για κάθε καταχώριση με το ίδιο ευρετήριο με df. Μεταβιβάζοντάς το στο df[] θα επιστρέψει αυτές τις καταχωρήσεις μόνο όταν οι αντίστοιχες τιμές τους είναι True.

Υπάρχουν 194 κλειδιά σε df[] για τα οποία df.source == 'GRCh38'.

Όπως είδαμε προηγουμένως, υπάρχουν επίσης 194 μοναδικές τιμές στο seqid στήλη, που σημαίνει κάθε καταχώριση σε gdf αντιστοιχεί σε ένα συγκεκριμένο seqid.

Στη συνέχεια επιλέγουμε τυχαία 10 καταχωρήσεις με το sample μέθοδος για πιο προσεκτική ματιά.

Μπορείτε να δείτε ότι οι μη συναρμολογημένες αλληλουχίες είναι τύπου supercontig ενώ οι άλλες είναι χρωμοσώματος. Για να υπολογίσουμε το κλάσμα του γονιδιώματος που είναι ελλιπές, πρέπει πρώτα να γνωρίζουμε το μήκος ολόκληρου του γονιδιώματος, που είναι το άθροισμα των μηκών όλων των σειρών.

In [90]: gdf = gdf.copy() In [91]: gdf['length'] = gdf.end - gdf.start + 1 In [93]: gdf.head() Out[93]: seqid source type start end score strand phase attributes length 0 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_00000... 248956421 235068 10 GRCh38 chromosome 1 133797422 . . . ID=chromosome:10;Alias=CM000672.2,chr10,NC_000... 133797421 328938 11 GRCh38 chromosome 1 135086622 . . . ID=chromosome:11;Alias=CM000673.2,chr11,NC_000... 135086621 483370 12 GRCh38 chromosome 1 133275309 . . . ID=chromosome:12;Alias=CM000674.2,chr12,NC_000... 133275308 634486 13 GRCh38 chromosome 1 114364328 . . . ID=chromosome:13;Alias=CM000675.2,chr13,NC_000... 114364327 In [97]: gdf.length.sum() Out[97]: 3096629532 In [99]: chrs = [str(_) for _ in range(1, 23)] + ['X', 'Y', 'MT'] In [101]: gdf[-gdf.seqid.isin(chrs)].length.sum() / gdf.length.sum() Out[101]: 0.0037021917421198327

Στο παραπάνω απόσπασμα, δημιουργήσαμε ένα αντίγραφο του gdf με .copy(). Διαφορετικά, το πρωτότυπο gdf είναι απλώς ένα κομμάτι df, και η άμεση τροποποίησή του θα είχε ως αποτέλεσμα SettingWithCopyWarning (βλέπω εδώ Για περισσότερες πληροφορίες).

Στη συνέχεια υπολογίζουμε τη διάρκεια κάθε καταχώρισης και την προσθέτουμε ξανά στο gdf ως μια νέα στήλη που ονομάζεται 'μήκος'. Το συνολικό μήκος είναι περίπου 3,1 δισεκατομμύρια και το κλάσμα των μη συναρμολογημένων ακολουθιών είναι περίπου 0,37%.

Δείτε πώς λειτουργεί ο τεμαχισμός στις δύο τελευταίες εντολές.

Πρώτον, δημιουργούμε μια λίστα με χορδές που καλύπτει όλες τις σειρές καλά συναρμολογημένων ακολουθιών, που είναι όλα τα χρωμοσώματα και τα μιτοχόνδρια. Στη συνέχεια χρησιμοποιούμε το isin μέθοδος για το φιλτράρισμα όλων των καταχωρήσεων των οποίων τα seqid βρίσκονται στο chrs λίστα.

Στην αρχή του ευρετηρίου προστίθεται ένα σύμβολο μείον (-) για να αντιστραφεί η επιλογή, γιατί πραγματικά θέλουμε ό, τι είναι δεν στη λίστα (δηλαδή θέλουμε τα μη συναρμολογημένα να ξεκινούν με KI και GL)…

Σημείωση: Δεδομένου ότι οι συναρμολογημένες και μη συναρμολογημένες ακολουθίες διακρίνονται από τη στήλη τύπου, η τελευταία γραμμή μπορεί εναλλακτικά να ξαναγραφεί ως εξής για να ληφθούν τα ίδια αποτελέσματα.

gdf[(gdf['type'] == 'supercontig')].length.sum() / gdf.length.sum()

Πόσα γονίδια υπάρχουν;

Εδώ επικεντρωνόμαστε στις καταχωρίσεις από το πηγαίο σύνολο, την Αβάνα και το ensembl_havana, καθώς βρίσκονται όπου ανήκουν οι περισσότερες καταχωρήσεις σχολιασμού.

In [109]: edf = df[df.source.isin(['ensembl', 'havana', 'ensembl_havana'])] In [111]: edf.sample(10) Out[111]: seqid source type start end score strand phase attributes 915996 16 havana CDS 27463541 27463592 . - 2 ID=CDS:ENSP00000457449;Parent=transcript:ENST0... 2531429 X havana exon 41196251 41196359 . + . Parent=transcript:ENST00000462850;Name=ENSE000... 1221944 19 ensembl_havana CDS 5641740 5641946 . + 0 ID=CDS:ENSP00000467423;Parent=transcript:ENST0... 243070 10 havana exon 13116267 13116340 . + . Parent=transcript:ENST00000378764;Name=ENSE000... 2413583 8 ensembl_havana exon 144359184 144359423 . + . Parent=transcript:ENST00000530047;Name=ENSE000... 2160496 6 havana exon 111322569 111322678 . - . Parent=transcript:ENST00000434009;Name=ENSE000... 839952 15 havana exon 76227713 76227897 . - . Parent=transcript:ENST00000565910;Name=ENSE000... 957782 16 ensembl_havana exon 67541653 67541782 . + . Parent=transcript:ENST00000379312;Name=ENSE000... 1632979 21 ensembl_havana exon 37840658 37840709 . - . Parent=transcript:ENST00000609713;Name=ENSE000... 1953399 4 havana exon 165464390 165464586 . + . Parent=transcript:ENST00000511992;Name=ENSE000... In [123]: edf.type.value_counts() Out[123]: exon 1180596 CDS 704604 five_prime_UTR 142387 three_prime_UTR 133938 transcript 96375 gene 42470 processed_transcript 28228 ... Name: type, dtype: int64

Το isin Η μέθοδος χρησιμοποιείται ξανά για φιλτράρισμα.

Στη συνέχεια, μια γρήγορη μέτρηση τιμών δείχνει ότι η πλειονότητα των καταχωρήσεων είναι εξόνιο, κωδικοποιητική ακολουθία (CDS) και μη μεταφρασμένη περιοχή (UTR).

Αυτά είναι στοιχεία υπο-γονιδίων, αλλά ψάχνουμε κυρίως για τον αριθμό των γονιδίων. Όπως φαίνεται, υπάρχουν 42.470, αλλά θέλουμε να μάθουμε περισσότερα.

Συγκεκριμένα, ποια είναι τα ονόματά τους και τι κάνουν; Για να απαντήσουμε σε αυτές τις ερωτήσεις, πρέπει να εξετάσουμε προσεκτικά τις πληροφορίες στη στήλη χαρακτηριστικών.

In [127]: ndf = edf[edf.type == 'gene'] In [173]: ndf = ndf.copy() In [133]: ndf.sample(10).attributes.values Out[133]: array(['ID=gene:ENSG00000228611;Name=HNF4GP1;biotype=processed_pseudogene;description=hepatocyte nuclear factor 4 gamma pseudogene 1 [Source:HGNC Symbol%3BAcc:HGNC:35417];gene_id=ENSG00000228611;havana_gene=OTTHUMG00000016986;havana_version=2;logic_name=havana;version=2', 'ID=gene:ENSG00000177189;Name=RPS6KA3;biotype=protein_coding;description=ribosomal protein S6 kinase A3 [Source:HGNC Symbol%3BAcc:HGNC:10432];gene_id=ENSG00000177189;havana_gene=OTTHUMG00000021231;havana_version=5;logic_name=ensembl_havana_gene;version=12', 'ID=gene:ENSG00000231748;Name=RP11-227H15.5;biotype=antisense;gene_id=ENSG00000231748;havana_gene=OTTHUMG00000018373;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000227426;Name=VN1R33P;biotype=unitary_pseudogene;description=vomeronasal 1 receptor 33 pseudogene [Source:HGNC Symbol%3BAcc:HGNC:37353];gene_id=ENSG00000227426;havana_gene=OTTHUMG00000154474;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000087250;Name=MT3;biotype=protein_coding;description=metallothionein 3 [Source:HGNC Symbol%3BAcc:HGNC:7408];gene_id=ENSG00000087250;havana_gene=OTTHUMG00000133282;havana_version=3;logic_name=ensembl_havana_gene;version=8', 'ID=gene:ENSG00000177108;Name=ZDHHC22;biotype=protein_coding;description=zinc finger DHHC-type containing 22 [Source:HGNC Symbol%3BAcc:HGNC:20106];gene_id=ENSG00000177108;havana_gene=OTTHUMG00000171575;havana_version=3;logic_name=ensembl_havana_gene;version=5', 'ID=gene:ENSG00000249784;Name=SCARNA22;biotype=scaRNA;description=small Cajal body-specific RNA 22 [Source:HGNC Symbol%3BAcc:HGNC:32580];gene_id=ENSG00000249784;logic_name=ncrna;version=1', 'ID=gene:ENSG00000079101;Name=CLUL1;biotype=protein_coding;description=clusterin like 1 [Source:HGNC Symbol%3BAcc:HGNC:2096];gene_id=ENSG00000079101;havana_gene=OTTHUMG00000178252;havana_version=7;logic_name=ensembl_havana_gene;version=16', 'ID=gene:ENSG00000229224;Name=AC105398.3;biotype=antisense;gene_id=ENSG00000229224;havana_gene=OTTHUMG00000152025;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000255552;Name=LY6G6E;biotype=protein_coding;description=lymphocyte antigen 6 complex%2C locus G6E (pseudogene) [Source:HGNC Symbol%3BAcc:HGNC:13934];gene_id=ENSG00000255552;havana_gene=OTTHUMG00000166419;havana_version=1;logic_name=ensembl_havana_gene;version=7'], dtype=object)

Διαμορφώνονται ως λίστα ζευγών τιμών-διαχωρισμένων με ερωτηματικά. Οι πληροφορίες που μας ενδιαφέρουν περισσότερο είναι το όνομα του γονιδίου, το αναγνωριστικό γονιδίου και η περιγραφή και θα τα εξαγάγουμε με κανονική έκφραση (regex).

import re RE_GENE_NAME = re.compile(r'Name=(?P.+?);') def extract_gene_name(attributes_str): res = RE_GENE_NAME.search(attributes_str) return res.group('gene_name') ndf['gene_name'] = ndf.attributes.apply(extract_gene_name)

Πρώτον, εξάγουμε τα ονόματα των γονιδίων.

Στο regex Name=(?P.+?); , +? χρησιμοποιείται αντί για + γιατί θέλουμε να μην είναι άπληστοι και να αφήσουμε την αναζήτηση να σταματήσει στο πρώτο ερωτηματικό. Διαφορετικά, το αποτέλεσμα θα ταιριάζει με το τελευταίο ερωτηματικό.

Επίσης, το regex καταρτίζεται για πρώτη φορά με re.compile αντί να χρησιμοποιείται απευθείας όπως στο re.search για καλύτερη απόδοση γιατί θα το εφαρμόσουμε σε χιλιάδες συμβολοσειρές χαρακτηριστικών.

extract_gene_name χρησιμεύει ως βοηθητική συνάρτηση που πρέπει να χρησιμοποιείται στο pd.apply, η οποία είναι η μέθοδος που χρησιμοποιείται όταν μια συνάρτηση πρέπει να εφαρμόζεται σε κάθε καταχώρηση ενός πλαισίου δεδομένων ή μιας σειράς.

Σε αυτήν τη συγκεκριμένη περίπτωση, θέλουμε να εξαγάγουμε το όνομα γονιδίου για κάθε καταχώριση στο ndf.attributes και να προσθέσουμε τα ονόματα πίσω στο ndf σε μια νέα στήλη που ονομάζεται gene_name.

Τα αναγνωριστικά γονιδίων και η περιγραφή εξάγονται με παρόμοιο τρόπο.

RE_GENE_ID = re.compile(r'gene_id=(?PENSG.+?);') def extract_gene_id(attributes_str): res = RE_GENE_ID.search(attributes_str) return res.group('gene_id') ndf['gene_id'] = ndf.attributes.apply(extract_gene_id) RE_DESC = re.compile('description=(?P.+?);') def extract_description(attributes_str): res = RE_DESC.search(attributes_str) if res is None: return '' else: return res.group('desc') ndf['desc'] = ndf.attributes.apply(extract_description)

Το regex για RE_GENE_ID είναι λίγο πιο συγκεκριμένη αφού γνωρίζουμε ότι κάθε gene_id πρέπει να ξεκινά με ENSG, όπου ENS που σημαίνει σύνολο και G σημαίνει γονίδιο.

Για καταχωρήσεις που δεν έχουν καμία περιγραφή, θα επιστρέψουμε μια κενή συμβολοσειρά. Μετά την εξαγωγή όλων, δεν θα χρησιμοποιούμε πλέον τη στήλη χαρακτηριστικών, οπότε ας το αφήσουμε για να διατηρήσουμε τα πράγματα ωραία και καθαρά με τη μέθοδο .drop:

In [224]: ndf.drop('attributes', axis=1, inplace=True) In [225]: ndf.head() Out[225]: seqid source type start end score strand phase gene_id gene_name desc 16 1 havana gene 11869 14409 . + . ENSG00000223972 DDX11L1 DEAD/H-box helicase 11 like 1 [Source:HGNC Sym... 28 1 havana gene 14404 29570 . - . ENSG00000227232 WASH7P WAS protein family homolog 7 pseudogene [Sourc... 71 1 havana gene 52473 53312 . + . ENSG00000268020 OR4G4P olfactory receptor family 4 subfamily G member... 74 1 havana gene 62948 63887 . + . ENSG00000240361 OR4G11P olfactory receptor family 4 subfamily G member... 77 1 ensembl_havana gene 69091 70008 . + . ENSG00000186092 OR4F5 olfactory receptor family 4 subfamily F member...

Στην παραπάνω κλήση, attributes υποδεικνύει τη συγκεκριμένη στήλη που θέλουμε να αφήσουμε.

axis=1 σημαίνει ότι ρίχνουμε μια στήλη αντί για μια σειρά (axis=0 από προεπιλογή).

inplace=True σημαίνει ότι η πτώση λειτουργεί στο ίδιο το DataFrame αντί να επιστρέφει ένα νέο αντίγραφο με καθορισμένη στήλη.

Ένα γρήγορο .head Η εμφάνιση δείχνει ότι η στήλη χαρακτηριστικών έχει πράγματι εξαφανιστεί και τρεις νέες στήλες: gene_name, gene_id και desc έχουν προστεθεί.

Από περιέργεια, ας δούμε αν όλα gene_id και gene_name είναι μοναδικά:

In [232]: ndf.shape Out[232]: (42470, 11) In [233]: ndf.gene_id.unique().shape Out[233]: (42470,) In [234]: ndf.gene_name.unique().shape Out[234]: (42387,)

Παραδόξως, ο αριθμός των γονιδιακών ονομάτων είναι μικρότερος από αυτόν των γονιδιακών αναγνωριστικών, υποδεικνύοντας ότι ορισμένα γονιδιακά ονόματα πρέπει να αντιστοιχούν σε πολλαπλά γονίδια ID. Ας μάθουμε τι είναι.

In [243]: count_df = ndf.groupby('gene_name').count().ix[:, 0].sort_values().ix[::-1] In [244]: count_df.head(10) Out[244]: gene_name SCARNA20 7 SCARNA16 6 SCARNA17 5 SCARNA15 4 SCARNA21 4 SCARNA11 4 Clostridiales-1 3 SCARNA4 3 C1QTNF9B-AS1 2 C11orf71 2 Name: seqid, dtype: int64 In [262]: count_df[count_df > 1].shape Out[262]: (63,) In [263]: count_df.shape Out[263]: (42387,) In [264]: count_df[count_df > 1].shape[0] / count_df.shape[0] Out[264]: 0.0014863047632528842

Ομαδοποιούμε όλες τις καταχωρίσεις με την τιμή gene_name και μετά μετράμε τον αριθμό των αντικειμένων σε κάθε ομάδα με .count().

Εάν ελέγξετε την έξοδο από ndf.groupby('gene_name').count(), όλες οι στήλες υπολογίζονται για κάθε ομάδα, αλλά οι περισσότερες από αυτές έχουν τις ίδιες τιμές.

Λάβετε υπόψη ότι οι τιμές NA δεν θα λαμβάνονται υπόψη κατά την καταμέτρηση, επομένως λάβετε μόνο το πλήθος της πρώτης στήλης, seqid (χρησιμοποιούμε .ix[:, 0] για να διασφαλίσουμε ότι δεν υπάρχουν τιμές NA).

Στη συνέχεια, ταξινομήστε τις τιμές μέτρησης με .sort_values και αντιστρέψτε τη σειρά με .ix[::-1].

Στο αποτέλεσμα, ένα όνομα γονιδίου μπορεί να κοινοποιηθεί σε έως και επτά αναγνωριστικά γονιδίων.

In [255]: ndf[ndf.gene_name == 'SCARNA20'] Out[255]: seqid source type start end score strand phase gene_id gene_name desc 179399 1 ensembl gene 171768070 171768175 . + . ENSG00000253060 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 201037 1 ensembl gene 204727991 204728106 . + . ENSG00000251861 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 349203 11 ensembl gene 8555016 8555146 . + . ENSG00000252778 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 718520 14 ensembl gene 63479272 63479413 . + . ENSG00000252800 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 837233 15 ensembl gene 75121536 75121666 . - . ENSG00000252722 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 1039874 17 ensembl gene 28018770 28018907 . + . ENSG00000251818 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 1108215 17 ensembl gene 60231516 60231646 . - . ENSG00000252577 SCARNA20 small Cajal body-specific RNA 20 [Source:HGNC Symbol%3BAcc:HGNC:32578]

Μια πιο προσεκτική ματιά σε όλα τα γονίδια SCARNA20 δείχνει ότι είναι πράγματι όλα διαφορετικά.

Ενώ έχουν το ίδιο όνομα, βρίσκονται σε διαφορετικές θέσεις του γονιδιώματος.

Ωστόσο, οι περιγραφές τους δεν φαίνονται πολύ χρήσιμες για τη διάκρισή τους.

Παράδειγμα πελάτη node js rest

Το σημείο εδώ είναι να γνωρίζουμε ότι τα ονόματα των γονιδίων δεν είναι μοναδικά για όλα τα αναγνωριστικά γονιδίων και περίπου 0,15% των ονομάτων που μοιράζονται πολλά γονίδια.

Πόσο καιρό είναι ένα τυπικό γονίδιο;

Παρόμοια με αυτό που κάναμε όταν προσπαθούσαμε να κατανοήσουμε την ατελή του γονιδιώματος, μπορούμε εύκολα να προσθέσουμε ένα length στήλη σε ndf:

In [277]: ndf['length'] = ndf.end - ndf.start + 1 In [278]: ndf.length.describe() Out[278]: count 4.247000e+04 mean 3.583348e+04 std 9.683485e+04 min 8.000000e+00 25% 8.840000e+02 50% 5.170500e+03 75% 3.055200e+04 max 2.304997e+06 Name: length, dtype: float64

.describe() υπολογίζει μερικά απλά στατιστικά στοιχεία με βάση τις τιμές μήκους:

• Το μέσο μήκος ενός γονιδίου είναι περίπου 36.000 βάσεις

• Το μέσο μήκος ενός γονιδίου έχει μήκος περίπου 5.200 βάσεις

• Τα ελάχιστα και μέγιστα μήκη γονιδίου έχουν μήκος περίπου οκτώ και 2,3 εκατομμύρια βάσεις, αντίστοιχα.

Επειδή ο μέσος όρος είναι πολύ μεγαλύτερος από τον διάμεσο, συνεπάγεται ότι η κατανομή μήκους είναι λοξή προς τα δεξιά. Για να έχουμε μια πιο συγκεκριμένη εμφάνιση, ας σχεδιάσουμε τη διανομή.

import matplotlib as plt ndf.length.plot(kind='hist', bins=50, logy=True) plt.show()

Το Pandas παρέχει μια απλή διεπαφή στο matplotlib για να κάνει τη σχεδίαση πολύ βολική με DataFrames ή σειρές.

Σε αυτήν την περίπτωση, λέει ότι θέλουμε μια γραφική παράσταση ιστογράμματος (kind='hist') με 50 κάδους και αφήστε τον άξονα y να βρίσκεται σε κλίμακα καταγραφής (logy=True).

Από το ιστόγραμμα, μπορούμε να δούμε ότι η πλειονότητα των γονιδίων βρίσκεται στον πρώτο κάδο. Ωστόσο, ορισμένα μήκη γονιδίων μπορεί να είναι πάνω από δύο εκατομμύρια βάσεις. Ας μάθουμε τι είναι:

In [39]: ndf[ndf.length > 2e6].sort_values('length').ix[::-1] Out[39]: seqid source type start end score strand phase gene_name gene_id desc length 2309345 7 ensembl_havana gene 146116002 148420998 . + . CNTNAP2 ENSG00000174469 contactin associated protein-like 2 [Source:HG... 2304997 2422510 9 ensembl_havana gene 8314246 10612723 . - . PTPRD ENSG00000153707 protein tyrosine phosphatase%2C receptor type ... 2298478 2527169 X ensembl_havana gene 31097677 33339441 . - . DMD ENSG00000198947 dystrophin [Source:HGNC Symbol%3BAcc:HGNC:2928] 2241765 440886 11 ensembl_havana gene 83455012 85627922 . - . DLG2 ENSG00000150672 discs large MAGUK scaffold protein 2 [Source:H... 2172911 2323457 8 ensembl_havana gene 2935353 4994972 . - . CSMD1 ENSG00000183117 CUB and Sushi multiple domains 1 [Source:HGNC ... 2059620 1569914 20 ensembl_havana gene 13995369 16053197 . + . MACROD2 ENSG00000172264 MACRO domain containing 2 [Source:HGNC Symbol%... 2057829

Όπως μπορείτε να δείτε, το μακρύτερο γονίδιο ονομάζεται CNTNAP2, το οποίο είναι συντομότερο για την πρωτεΐνη που σχετίζεται με την πρωτεΐνη 2. Σύμφωνα με το σελίδα wikipedia ,

Αυτό το γονίδιο περιλαμβάνει σχεδόν το 1,6% του χρωμοσώματος 7 και είναι ένα από τα μεγαλύτερα γονίδια στο ανθρώπινο γονιδίωμα.

Πράγματι! Μόλις το επιβεβαιώσαμε. Αντίθετα, τι γίνεται με τα μικρότερα γονίδια; Αποδεικνύεται ότι μπορούν να είναι τόσο μικρές όσο οκτώ βάσεις.

In [40]: ndf.sort_values('length').head() Out[40]: seqid source type start end score strand phase gene_name gene_id desc length 682278 14 havana gene 22438547 22438554 . + . TRDD1 ENSG00000223997 T cell receptor delta diversity 1 [Source:HGNC... 8 682282 14 havana gene 22439007 22439015 . + . TRDD2 ENSG00000237235 T cell receptor delta diversity 2 [Source:HGNC... 9 2306836 7 havana gene 142786213 142786224 . + . TRBD1 ENSG00000282431 T cell receptor beta diversity 1 [Source:HGNC ... 12 682286 14 havana gene 22449113 22449125 . + . TRDD3 ENSG00000228985 T cell receptor delta diversity 3 [Source:HGNC... 13 1879625 4 havana gene 10238213 10238235 . - . AC006499.9 ENSG00000271544 23

Το μήκος των δύο ακραίων περιπτώσεων είναι πέντε τάξεις μεγέθους (2,3 εκατομμύρια έναντι 8), το οποίο είναι τεράστιο και μπορεί να είναι ένδειξη του επιπέδου ποικιλομορφίας της ζωής.

Ένα μεμονωμένο γονίδιο μπορεί να μεταφραστεί σε πολλές διαφορετικές πρωτεΐνες μέσω μιας διαδικασίας που ονομάζεται εναλλακτική σύνδεση, κάτι που δεν έχουμε διερευνήσει. Τέτοιες πληροφορίες βρίσκονται επίσης εντός του αρχείου GFF3, αλλά εκτός του πεδίου εφαρμογής αυτής της ανάρτησης.

Κατανομή γονιδίων μεταξύ χρωμοσωμάτων

Το τελευταίο πράγμα που θα ήθελα να συζητήσω είναι η κατανομή γονιδίων μεταξύ των χρωμοσωμάτων, η οποία χρησιμεύει επίσης ως παράδειγμα για την εισαγωγή του .merge μέθοδος συνδυασμού δύο DataFrames. Διαισθητικά, μακρύτερα χρωμοσώματα πιθανώς φιλοξενούν περισσότερα γονίδια. Ας δούμε αν αυτό είναι αλήθεια.

In [53]: ndf = ndf[ndf.seqid.isin(chrs)] In [54]: chr_gene_counts = ndf.groupby('seqid').count().ix[:, 0].sort_values().ix[::-1] Out[54]: chr_gene_counts seqid 1 3902 2 2806 11 2561 19 2412 17 2280 3 2204 6 2154 12 2140 7 2106 5 2002 16 1881 X 1852 4 1751 9 1659 8 1628 10 1600 15 1476 14 1449 22 996 20 965 13 872 18 766 21 541 Y 436 Name: source, dtype: int64

Δανειστήκαμε το chrs μεταβλητή από την προηγούμενη ενότητα και το χρησιμοποίησε για να φιλτράρει τις μη συναρμολογημένες ακολουθίες. Με βάση την έξοδο, το μεγαλύτερο χρωμόσωμα 1 έχει όντως τα περισσότερα γονίδια. Ενώ το χρωμόσωμα Υ έχει τον μικρότερο αριθμό γονιδίων, δεν είναι το μικρότερο χρωμόσωμα.

Σημειώστε ότι δεν φαίνεται να υπάρχουν γονίδια στο μιτοχόνδριο (MT), κάτι που δεν ισχύει.

Λίγο πιο φιλτράρισμα στο πρώτο DataFrame df επέστρεψε από pd.read_csv δείχνει ότι όλα τα γονίδια MT προέρχονται από την πηγή insdc (τα οποία είχαν φιλτραριστεί πριν κατά τη δημιουργία edf όπου εξετάσαμε μόνο πηγές αβάνας, ensembl ή ensembl_havana).

In [134]: df[(df.type == 'gene') & (df.seqid == 'MT')] Out[134]: seqid source type start end score strand phase attributes 2514003 MT insdc gene 648 1601 . + . ID=gene:ENSG00000211459;Name=MT-RNR1;biotype=M... 2514009 MT insdc gene 1671 3229 . + . ID=gene:ENSG00000210082;Name=MT-RNR2;biotype=M... 2514016 MT insdc gene 3307 4262 . + . ID=gene:ENSG00000198888;Name=MT-ND1;biotype=pr... 2514029 MT insdc gene 4470 5511 . + . ID=gene:ENSG00000198763;Name=MT-ND2;biotype=pr... 2514048 MT insdc gene 5904 7445 . + . ID=gene:ENSG00000198804;Name=MT-CO1;biotype=pr... 2514058 MT insdc gene 7586 8269 . + . ID=gene:ENSG00000198712;Name=MT-CO2;biotype=pr... 2514065 MT insdc gene 8366 8572 . + . ID=gene:ENSG00000228253;Name=MT-ATP8;biotype=p... 2514069 MT insdc gene 8527 9207 . + . ID=gene:ENSG00000198899;Name=MT-ATP6;biotype=p... 2514073 MT insdc gene 9207 9990 . + . ID=gene:ENSG00000198938;Name=MT-CO3;biotype=pr... 2514080 MT insdc gene 10059 10404 . + . ID=gene:ENSG00000198840;Name=MT-ND3;biotype=pr... 2514087 MT insdc gene 10470 10766 . + . ID=gene:ENSG00000212907;Name=MT-ND4L;biotype=p... 2514091 MT insdc gene 10760 12137 . + . ID=gene:ENSG00000198886;Name=MT-ND4;biotype=pr... 2514104 MT insdc gene 12337 14148 . + . ID=gene:ENSG00000198786;Name=MT-ND5;biotype=pr... 2514108 MT insdc gene 14149 14673 . - . ID=gene:ENSG00000198695;Name=MT-ND6;biotype=pr... 2514115 MT insdc gene 14747 15887 . + . ID=gene:ENSG00000198727;Name=MT-CYB;biotype=pr...

Αυτό το παράδειγμα δείχνει επίσης πώς να συνδυάσετε δύο συνθήκες κατά τη διάρκεια του φιλτραρίσματος με &; ο λογικός τελεστής για 'ή' θα είναι |.

Σημειώστε ότι απαιτούνται παρενθέσεις γύρω από κάθε συνθήκη και αυτό το μέρος της σύνταξης στο Pandas είναι διαφορετικό από το Python, το οποίο θα ήταν κυριολεκτικό and και or.

Στη συνέχεια, ας δανειστεί το gdf DataFrame από την προηγούμενη ενότητα ως πηγή για το μήκος κάθε χρωμοσώματος:

In [61]: gdf = gdf[gdf.seqid.isin(chrs)] In [62]: gdf.drop(['start', 'end', 'score', 'strand', 'phase' ,'attributes'], axis=1, inplace=True) In [63]: gdf.sort_values('length').ix[::-1] Out[63]: seqid source type length 0 1 GRCh38 chromosome 248956422 1364641 2 GRCh38 chromosome 242193529 1705855 3 GRCh38 chromosome 198295559 1864567 4 GRCh38 chromosome 190214555 1964921 5 GRCh38 chromosome 181538259 2080148 6 GRCh38 chromosome 170805979 2196981 7 GRCh38 chromosome 159345973 2514125 X GRCh38 chromosome 156040895 2321361 8 GRCh38 chromosome 145138636 2416560 9 GRCh38 chromosome 138394717 328938 11 GRCh38 chromosome 135086622 235068 10 GRCh38 chromosome 133797422 483370 12 GRCh38 chromosome 133275309 634486 13 GRCh38 chromosome 114364328 674767 14 GRCh38 chromosome 107043718 767312 15 GRCh38 chromosome 101991189 865053 16 GRCh38 chromosome 90338345 990810 17 GRCh38 chromosome 83257441 1155977 18 GRCh38 chromosome 80373285 1559144 20 GRCh38 chromosome 64444167 1201561 19 GRCh38 chromosome 58617616 2594560 Y GRCh38 chromosome 54106423 1647482 22 GRCh38 chromosome 50818468 1616710 21 GRCh38 chromosome 46709983 2513999 MT GRCh38 chromosome 16569

Οι στήλες που δεν σχετίζονται με την ανάλυση απορρίπτονται για λόγους σαφήνειας.

Ναι, .drop μπορεί επίσης να λάβει μια λίστα στηλών και να τις ρίξει εντελώς σε μία λειτουργία.

Σημειώστε ότι η σειρά με seqid MT είναι ακόμα εκεί. θα το επιστρέψουμε αργότερα. Η επόμενη λειτουργία που θα εκτελέσουμε είναι η συγχώνευση των δύο συνόλων δεδομένων με βάση τις τιμές του seqid.

In [73]: cdf = chr_gene_counts.to_frame(name='gene_count').reset_index() In [75]: cdf.head(2) Out[75]: seqid gene_count 0 1 3902 1 2 2806 In [78]: merged = gdf.merge(cdf, on='seqid') In [79]: merged Out[79]: seqid source type length gene_count 0 1 GRCh38 chromosome 248956422 3902 1 10 GRCh38 chromosome 133797422 1600 2 11 GRCh38 chromosome 135086622 2561 3 12 GRCh38 chromosome 133275309 2140 4 13 GRCh38 chromosome 114364328 872 5 14 GRCh38 chromosome 107043718 1449 6 15 GRCh38 chromosome 101991189 1476 7 16 GRCh38 chromosome 90338345 1881 8 17 GRCh38 chromosome 83257441 2280 9 18 GRCh38 chromosome 80373285 766 10 19 GRCh38 chromosome 58617616 2412 11 2 GRCh38 chromosome 242193529 2806 12 20 GRCh38 chromosome 64444167 965 13 21 GRCh38 chromosome 46709983 541 14 22 GRCh38 chromosome 50818468 996 15 3 GRCh38 chromosome 198295559 2204 16 4 GRCh38 chromosome 190214555 1751 17 5 GRCh38 chromosome 181538259 2002 18 6 GRCh38 chromosome 170805979 2154 19 7 GRCh38 chromosome 159345973 2106 20 8 GRCh38 chromosome 145138636 1628 21 9 GRCh38 chromosome 138394717 1659 22 X GRCh38 chromosome 156040895 1852 23 Y GRCh38 chromosome 54106423 436

Από chr_gene_counts εξακολουθεί να είναι ένα αντικείμενο Σειράς, το οποίο δεν υποστηρίζει μια λειτουργία συγχώνευσης, πρέπει πρώτα να το μετατρέψουμε σε αντικείμενο DataFrame με .to_frame.

.reset_index() μετατρέπει το αρχικό ευρετήριο (δηλαδή seqid) σε μια νέα στήλη και επαναφέρει το τρέχον ευρετήριο ως αυξητικούς αριθμούς με βάση 0.

Η έξοδος από cdf.head(2) δείχνει πώς μοιάζει. Στη συνέχεια, χρησιμοποιήσαμε το .merge μέθοδος συνδυασμού των δύο DataFrame στη στήλη seqid (on='seqid').

Μετά τη συγχώνευση gdf και cdf, η MT λείπει η καταχώριση. Αυτό συμβαίνει επειδή, από προεπιλογή, .merge λειτουργεί μια εσωτερική ένωση, ενώ η αριστερή ένωση, η δεξιά ένωση ή η εξωτερική ένωση είναι διαθέσιμες συντονίζοντας το how παράμετρος.

Ανατρέξτε στο τεκμηρίωση Για περισσότερες πληροφορίες.

Αργότερα, μπορεί να διαπιστώσετε ότι υπάρχει επίσης μια σχετική .join μέθοδος. .merge και .join είναι παρόμοια αλλά έχουν διαφορετικά API.

Σύμφωνα με τον αξιωματούχο τεκμηρίωση λέει

Η σχετική μέθοδος DataFrame.join, χρησιμοποιεί συγχώνευση εσωτερικά για το ευρετήριο-ευρετήριο και το ευρετήριο-στη στήλη (ες) ενώνει, αλλά συνδέεται στα ευρετήρια από προεπιλογή αντί να προσπαθεί να ενωθεί σε κοινές στήλες (η προεπιλεγμένη συμπεριφορά για συγχώνευση). Εάν συμμετέχετε στο ευρετήριο, ίσως θελήσετε να χρησιμοποιήσετε το DataFrame.join για να σώσετε τον εαυτό σας κάποια πληκτρολόγηση.

Βασικά, .merge είναι πιο γενικής χρήσης και χρησιμοποιείται από .join.

Τέλος, είμαστε έτοιμοι να υπολογίσουμε τη συσχέτιση μεταξύ του χρωμοσώματος length και gene_count.

In [81]: merged[['length', 'gene_count']].corr() Out[81]: length gene_count length 1.000000 0.728221 gene_count 0.728221 1.000000

Από προεπιλογή .corr υπολογίζει το Συσχέτιση Pearson μεταξύ όλων των ζευγών στηλών σε ένα πλαίσιο δεδομένων.

Στην περίπτωση αυτή έχουμε μόνο ένα ζεύγος στηλών και η συσχέτιση αποδεικνύεται θετική - 0,73.

Με άλλα λόγια, όσο μεγαλύτερο είναι το χρωμόσωμα, τόσο πιθανότερο είναι να έχουμε περισσότερα γονίδια. Ας σχεδιάσουμε επίσης τις δύο στήλες αφού ταξινομήσουμε τα ζεύγη τιμών κατά length:

ax = merged[['length', 'gene_count']].sort_values('length').plot(x='length', y='gene_count', style='o-') # add some margin to both ends of x axis xlim = ax.get_xlim() margin = xlim[0] * 0.1 ax.set_xlim([xlim[0] - margin, xlim[1] + margin]) # Label each point on the graph for (s, x, y) in merged[['seqid', 'length', 'gene_count']].sort_values('length').values: ax.text(x, y - 100, str(s))

Όπως φαίνεται στην παραπάνω εικόνα, παρόλο που είναι συνολικά θετικός συσχετισμός, δεν ισχύει για όλα τα χρωμοσώματα. Συγκεκριμένα, για το χρωμόσωμα 17, 16, 15, 14, 13, ο συσχετισμός είναι στην πραγματικότητα αρνητικός, πράγμα που σημαίνει ότι ο αριθμός των γονιδίων στο χρωμόσωμα μειώνεται καθώς αυξάνεται το μέγεθος του χρωμοσώματος.

Ευρήματα και μελλοντική έρευνα

Αυτό τελειώνει το σεμινάριό μας για το χειρισμό ενός αρχείου σχολιασμού για το ανθρώπινο γονιδίωμα σε μορφή GFF3 με τη στοίβα SciPy. Τα εργαλεία που χρησιμοποιήσαμε κυρίως περιλαμβάνουν IPython, Pandas και matplotlib. Κατά τη διάρκεια του σεμιναρίου, όχι μόνο μάθαμε μερικές από τις πιο κοινές και χρήσιμες λειτουργίες στο Pandas, αλλά και απαντήσαμε σε μερικές πολύ ενδιαφέρουσες ερωτήσεις σχετικά με το γονιδίωμά μας. Συνοψίζοντας:

1. Περίπου το 0,37% του ανθρώπινου γονιδιώματος είναι ακόμη ατελές, παρόλο που το πρώτο σχέδιο κυκλοφόρησε πριν από 15 χρόνια.
2. Υπάρχουν περίπου 42.000 γονίδια στο ανθρώπινο γονιδίωμα με βάση αυτό το συγκεκριμένο Αρχείο GFF3 συνηθίζαμε.
3. Το μήκος ενός γονιδίου μπορεί να κυμαίνεται από μερικές δεκάδες έως πάνω από δύο εκατομμύρια βάσεις.
4. Τα γονίδια δεν κατανέμονται ομοιόμορφα μεταξύ των χρωμοσωμάτων. Συνολικά, όσο μεγαλύτερο είναι το χρωμόσωμα, τόσο περισσότερα γονίδια φιλοξενεί, αλλά για ένα υποσύνολο των χρωμοσωμάτων, ο συσχετισμός μπορεί να είναι αρνητικός.

Το αρχείο GFF3 είναι πολύ πλούσιο σε πληροφορίες σχολιασμού και έχουμε μόλις γρατσουνίσει την επιφάνεια. Εάν ενδιαφέρεστε για περαιτέρω εξερεύνηση, ακολουθούν μερικές ερωτήσεις με τις οποίες μπορείτε να παίξετε:

1. Πόσες μεταγραφές έχει συνήθως ένα γονίδιο; Ποιο ποσοστό γονιδίων έχουν περισσότερα από 1 μεταγραφή;
2. Πόσες ισομορφές έχει συνήθως ένα αντίγραφο;
3. Πόσα εξόνια, CDS και UTR έχει συνήθως μια μεταγραφή; Τι μεγέθη είναι;
4. Είναι δυνατόν να κατηγοριοποιηθούν τα γονίδια με βάση τη λειτουργία τους όπως περιγράφεται στη στήλη περιγραφής;