Eine umfassende Einführung in Ihr Genom mit dem SciPy Stack

Die Bioinformatik ist ein interdisziplinäres Gebiet, das Methoden und Softwaretools zur Analyse und zum Verständnis biologischer Daten entwickelt.

Einfacher ausgedrückt kann man es sich einfach als Datenwissenschaft für die Biologie vorstellen.

Unter den vielen Arten biologischer Daten gehören Genomdaten zu den am häufigsten analysierten. Insbesondere mit der schnellen Weiterentwicklung der DNA-Sequenzierungstechnologien der nächsten Generation (NGS) ist das Volumen der Genomdaten exponentiell gewachsen. Laut Stephens, Zachary D et al. liegt die genomische Datenerfassung im Bereich von Exabyte pro Jahr.

In diesem Beitrag zeige ich ein Beispiel für die Analyse einer GFF3-Datei für das menschliche Genom mit dem SciPy Stack. Generic Feature Format Version 3 (GFF3) ist das aktuelle Standard-Textdateiformat zum Speichern genomischer Merkmale. Insbesondere erfahren Sie in diesem Beitrag, wie Sie den SciPy-Stack verwenden, um die folgenden Fragen zum menschlichen Genom zu beantworten:

1. Wie viel des Genoms ist unvollständig?
2. Wie viele Gene enthält das Genom?
3. Wie lang ist ein typisches Gen?
4. Wie sieht die Genverteilung zwischen den Chromosomen aus?

Die neueste GFF3-Datei für das menschliche Genom kann hier heruntergeladen werden. Die README-Datei, die sich in diesem Verzeichnis befindet, enthält eine kurze Beschreibung dieses Datenformats, und eine ausführlichere Spezifikation finden Sie hier.

Wir werden Pandas verwenden, eine Hauptkomponente des SciPy-Stacks, die schnelle, flexible und ausdrucksstarke Datenstrukturen bereitstellt, um die GFF3-Datei zu manipulieren und zu verstehen.

Aufstellen

Lassen Sie uns zuerst eine virtuelle Umgebung mit installiertem SciPy-Stack einrichten. Dieser Prozess kann zeitaufwändig sein, wenn er manuell aus der Quelle erstellt wird, da der Stack viele Pakete umfasst – von denen einige von externem FORTRAN- oder C-Code abhängen. Hier empfehle ich die Verwendung von Miniconda, was die Einrichtung sehr einfach macht.

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

Das Flag -b in der Bash-Zeile weist es an, im Batch-Modus ausgeführt zu werden. Nachdem die obigen Befehle verwendet wurden, um Miniconda erfolgreich zu installieren, starten Sie eine neue virtuelle Umgebung für die Genomik und installieren Sie dann den SciPy-Stack.

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

Beachten Sie, dass wir in diesem Beitrag nur die 3 Pakete angegeben haben, die wir verwenden werden. Wenn Sie möchten, dass alle Pakete im SciPy-Stack aufgelistet sind, hängen Sie sie einfach an das Ende des conda create Befehls an.

Wenn Sie den genauen Namen eines Pakets nicht kennen, versuchen Sie es mit conda search . Lassen Sie uns die virtuelle Umgebung aktivieren und IPython starten.

 source activate venv/ ipython

IPython ist ein deutlich leistungsstärkerer Ersatz für die standardmäßige Python-Interpreterschnittstelle, sodass alles, was Sie früher im standardmäßigen Python-Interpreter getan haben, auch in IPython ausgeführt werden kann. Ich empfehle jedem Python-Programmierer, der IPython noch nicht verwendet hat, es einmal auszuprobieren.

Laden Sie die Anmerkungsdatei herunter

Nachdem unser Setup nun abgeschlossen ist, laden wir die Annotationsdatei des menschlichen Genoms im GFF3-Format herunter.

Es ist etwa 37 MB groß, eine sehr kleine Datei im Vergleich zum Informationsgehalt eines menschlichen Genoms, das im Klartext etwa 3 GB groß ist. Das liegt daran, dass die GFF3-Datei nur die Annotation der Sequenzen enthält, während die Sequenzdaten normalerweise in einem anderen Dateiformat namens FASTA gespeichert werden. Wenn Sie interessiert sind, können Sie FASTA hier herunterladen, aber wir werden die Sequenzdaten in diesem Tutorial nicht verwenden.

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

Das vorangestellte ! teilt IPython mit, dass dies ein Shell-Befehl anstelle eines Python-Befehls ist. Allerdings kann IPython einige häufig verwendete Shell-Befehle wie ls , pwd , rm , mkdir , rmdir auch ohne vorangestelltes ! .

Wenn Sie sich den Kopf der GFF-Datei ansehen, werden Sie viele Zeilen mit Metadaten/Pragmas/Anweisungen sehen, die mit ## oder #! .

Laut README bedeutet ## , dass die Metadaten stabil sind, während #! bedeutet, es ist experimentell.

Später sehen Sie auch ### , eine weitere Direktive mit noch subtilerer Bedeutung, basierend auf der Spezifikation.

Für Menschen lesbare Kommentare sollten nach einem einzelnen # stehen. Der Einfachheit halber behandeln wir alle Zeilen, die mit # beginnen, als Kommentare und ignorieren sie einfach während unserer Analyse.

 ##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

Die erste Zeile gibt an, dass die in dieser Datei verwendete Version des GFF-Formats 3 ist.

Darauf folgen Zusammenfassungen aller Sequenzregionen. Wie wir später sehen werden, können solche Informationen auch im Body-Teil der Datei gefunden werden.

Die Zeilen beginnend mit #! zeigt Informationen über den bestimmten Build des Genoms, GRCh38.p7, auf das sich diese Anmerkungsdatei bezieht.

Das Genome Reference Consortium (GCR) ist ein internationales Konsortium, das Aktualisierungen und Verbesserungen mehrerer Referenzgenombaugruppen überwacht, einschließlich derjenigen für Menschen, Mäuse, Zebrafische und Hühner.

Beim Scannen dieser Datei sind hier die ersten Anmerkungszeilen.

 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 ...

Die Spalten sind seqid , Quelle , Typ , Start , Ende , Punktzahl , Strang , Phase , Attribute . Einige von ihnen sind sehr einfach zu verstehen. Nehmen wir als Beispiel die erste Zeile:

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

Dies ist die Annotation des ersten Chromosoms mit einer Seqid von 1, die von der ersten Base bis zur 24.895.622. Base beginnt.

Mit anderen Worten, das erste Chromosom ist etwa 25 Millionen Basen lang.

Unsere Analyse benötigt keine Informationen aus den drei Spalten mit einem Wert von . (dh Partitur, Strang und Phase), sodass wir sie vorerst einfach ignorieren können.

Die letzte Attributspalte besagt, dass Chromosom 1 auch drei Aliasnamen hat, nämlich CM000663.2, chr1 und NC_000001.11. So sieht eine GFF3-Datei im Grunde aus, aber wir werden sie nicht Zeile für Zeile untersuchen, also ist es an der Zeit, die gesamte Datei in Pandas zu laden.

Pandas eignet sich gut für den Umgang mit dem GFF3-Format, da es sich um eine tabulatorgetrennte Datei handelt, und Pandas bietet eine sehr gute Unterstützung für das Lesen von CSV-ähnlichen Dateien.

Beachten Sie, dass eine Ausnahme vom tabulatorgetrennten Format auftritt, wenn GFF3 ## ##FASTA enthält.

Gemäß der Spezifikation gibt ##FASTA das Ende eines Anmerkungsabschnitts an, dem eine oder mehrere Sequenzen im FASTA-Format (einem nicht durch Tabulatoren getrennten Format) folgen. Dies ist jedoch bei der GFF3-Datei, die wir analysieren werden, nicht der Fall.

 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='\t', comment='#', low_memory=False, header=None, names=col_names)

Die letzte Zeile oben lädt die gesamte GFF3-Datei mit der pandas.read_csv Methode.

Da es sich nicht um eine Standard-CSV-Datei handelt, müssen wir den Aufruf etwas anpassen.

Zuerst informieren wir Pandas über die Nichtverfügbarkeit von Header-Informationen in GFF3 mit header=None und geben dann den genauen Namen für jede Spalte mit names=col_names .

Wenn das Argument names nicht angegeben ist, verwendet Pandas inkrementelle Zahlen als Namen für jede Spalte.

sep='\t' teilt Pandas mit, dass die Spalten tabulatorgetrennt statt kommagetrennt sind. Der Wert von sep= kann tatsächlich ein regulärer Ausdruck (regex) sein. Dies wird praktisch, wenn die vorliegende Datei unterschiedliche Trennzeichen für jede Spalte verwendet (oh ja, das passiert). comment='#' bedeutet, dass Zeilen, die mit # beginnen, als Kommentare betrachtet und ignoriert werden.

compression='gzip' teilt Pandas mit, dass die Eingabedatei gzip-komprimiert ist.

Darüber hinaus verfügt pandas.read_csv über einen umfangreichen Satz von Parametern, mit denen verschiedene Arten von CSV-ähnlichen Dateiformaten gelesen werden können.

Der Typ des zurückgegebenen Werts ist ein DataFrame , die wichtigste Datenstruktur in Pandas, die zur Darstellung von 2D-Daten verwendet wird.

Pandas hat auch eine Series und Panel -Datenstruktur für 1D- bzw. 3D-Daten. Eine Einführung in die Datenstrukturen von Pandas finden Sie in der Dokumentation.

Werfen wir einen Blick auf die ersten paar Einträge mit der .head Methode.

 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

Die Ausgabe ist schön in einem tabellarischen Format formatiert, wobei längere Zeichenfolgen in der Attribute-Spalte teilweise durch ... ersetzt werden.

Mit pd.set_option('display.max_colwidth', -1) können Sie Pandas so einstellen, dass lange Zeichenfolgen nicht weggelassen werden. Darüber hinaus bietet Pandas viele Optionen, die angepasst werden können.

Lassen Sie uns als Nächstes einige grundlegende Informationen zu diesem Datenrahmen mit der .info -Methode abrufen.

 In [20]: df.info() <class 'pandas.core.frame.DataFrame'> 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

Dies zeigt, dass das GFF3 2.601.848 kommentierte Zeilen hat und jede Zeile neun Spalten hat.

Für jede Spalte werden auch ihre Datentypen angezeigt.

Dieser Start und dieses Ende sind vom Typ int64, ganze Zahlen, die Positionen im Genom darstellen.

Die anderen Spalten sind alle vom Typ object , was wahrscheinlich bedeutet, dass ihre Werte aus einer Mischung von Ganzzahlen, Floats und Strings bestehen.

Die Größe aller Informationen beträgt etwa 178,7+ MB im Speicher. Dies erweist sich als kompakter als die unkomprimierte Datei, die etwa 402 MB groß sein wird. Eine schnelle Überprüfung wird unten gezeigt.

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

Aus einer übergeordneten Sicht haben wir die gesamte GFF3-Datei in ein DataFrame-Objekt in Python geladen, und alle unsere folgenden Analysen basieren auf diesem einzelnen Objekt.

Sehen wir uns nun an, worum es bei der ersten Spalte seqid geht.

 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 ist eine Möglichkeit, auf Spaltendaten aus einem Datenrahmen zuzugreifen. Ein anderer Weg ist df['seqid'] , was eine allgemeinere Syntax ist, denn wenn der Spaltenname ein für Python reserviertes Schlüsselwort ist (z. B. class ) oder eine . oder Leerzeichen, der erste Weg ( df.seqid ) funktioniert nicht.

Die Ausgabe zeigt, dass es 194 eindeutige Seqids gibt, darunter Chromosomen 1 bis 22, X-, Y- und Mitochondrien (MT)-DNA sowie 169 andere Seqids.

Die mit KI und GL beginnenden Sequenzen sind DNA-Sequenzen – oder Gerüste – im Genom, die nicht erfolgreich in die Chromosomen eingebaut wurden.

Für diejenigen, die mit Genomik nicht vertraut sind, ist dies wichtig.

Obwohl der erste Entwurf des menschlichen Genoms vor mehr als 15 Jahren herauskam, ist das aktuelle menschliche Genom immer noch unvollständig. Die Schwierigkeit beim Zusammensetzen dieser Sequenzen ist größtenteils auf komplexe repetitive Regionen im Genom zurückzuführen.

Als nächstes werfen wir einen Blick auf die Quellspalte.

Die README besagt, dass die Quelle ein Freitextqualifizierer ist, der den Algorithmus oder die Betriebsprozedur beschreiben soll, die diese Funktion generiert hat.

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

Dies ist ein Beispiel für die Verwendung der value_counts Methode, die für eine schnelle Zählung kategorialer Variablen äußerst nützlich ist.

Aus dem Ergebnis können wir ersehen, dass es sieben mögliche Werte für diese Spalte gibt und die meisten Einträge in der GFF3-Datei von havana, ensembl und ensembl_havana stammen.

In diesem Beitrag erfahren Sie mehr über die Bedeutung dieser Quellen und die Beziehungen zwischen ihnen.

Um die Dinge einfach zu halten, konzentrieren wir uns auf Einträge aus den Quellen GRCh38, havana, ensembl und ensembl_havan.a.

Wie viel des Genoms ist unvollständig?

Die Informationen über jedes gesamte Chromosom befinden sich in den Einträgen der Quelle GRCh38, also lassen Sie uns zuerst den Rest herausfiltern und das gefilterte Ergebnis einer neuen Variablen 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...

Das Filtern ist in Pandas einfach.

Wenn Sie den aus dem Ausdruck df.source == 'GRCh38' ausgewerteten Wert untersuchen, handelt es sich um eine Reihe von True und False -Werten für jeden Eintrag mit demselben Index wie df . Wenn Sie es an df[] , werden nur die Einträge zurückgegeben, deren entsprechende Werte True sind.

Es gibt 194 Schlüssel in df[] , für die df.source == 'GRCh38' .

Wie wir zuvor gesehen haben, gibt es auch 194 eindeutige Werte in der seqid Spalte, was bedeutet, dass jeder Eintrag in gdf einer bestimmten seqid entspricht.

Dann wählen wir zufällig 10 Einträge mit der sample aus, um sie uns genauer anzusehen.

Sie können sehen, dass die nicht zusammengesetzten Sequenzen vom Typ Supercontig sind, während die anderen vom Chromosom sind. Um den Anteil des unvollständigen Genoms zu berechnen, müssen wir zunächst die Länge des gesamten Genoms kennen, die die Summe der Längen aller Seqids ist.

 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

Im obigen Snippet haben wir zuerst eine Kopie von gdf mit .copy() . Andernfalls ist das ursprüngliche gdf nur ein Teil von df , und eine direkte Änderung würde zu SettingWithCopyWarning (weitere Einzelheiten finden Sie hier).

Wir berechnen dann die Länge jedes Eintrags und fügen sie als neue Spalte mit dem Namen „Länge“ wieder zu gdf hinzu. Die Gesamtlänge beträgt etwa 3,1 Milliarden, und der Anteil nicht zusammengesetzter Sequenzen beträgt etwa 0,37 %.

So funktioniert das Slicing in den letzten beiden Befehlen.

Zuerst erstellen wir eine Liste von Zeichenfolgen, die alle Seqids gut zusammengesetzter Sequenzen abdeckt, die alle Chromosomen und Mitochondrien sind. Anschließend filtern wir mit der isin Methode alle Einträge, deren seqid in der chrs Liste stehen.

Ein Minuszeichen ( - ) wird am Anfang des Index hinzugefügt, um die Auswahl umzukehren, weil wir eigentlich alles wollen, was nicht in der Liste ist (dh wir wollen die unassemblierten, beginnend mit KI und GL)…

Hinweis: Da die zusammengesetzten und nicht zusammengesetzten Sequenzen durch die Typspalte unterschieden werden, kann die letzte Zeile alternativ wie folgt umgeschrieben werden, um die gleichen Ergebnisse zu erhalten.

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

Wie viele Gene gibt es?

Hier konzentrieren wir uns auf die Einträge aus den Quellen ensembl, havana und ensembl_havana, da dort die meisten Anmerkungseinträge hingehören.

 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

Zur Filterung wird wieder die isin Methode verwendet.

Dann zeigt eine schnelle Wertezählung, dass die Mehrzahl der Einträge Exon, kodierende Sequenz (CDS) und untranslatierte Region (UTR) sind.

Dies sind Subgenelemente, aber wir suchen hauptsächlich nach der Genzahl. Wie gezeigt, gibt es 42.470, aber wir wollen mehr wissen.

Genauer gesagt, wie heißen sie und was machen sie? Um diese Fragen zu beantworten, müssen wir uns die Informationen in der Attribute-Spalte genau ansehen.

 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)

Sie sind als durch Semikolons getrennte Liste von Tag-Wert-Paaren formatiert. Die Informationen, an denen wir am meisten interessiert sind, sind Genname, Gen-ID und Beschreibung, und wir werden sie mit regulären Ausdrücken (Regex) extrahieren.

 import re RE_GENE_NAME = re.compile(r'Name=(?P<gene_name>.+?);') 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)

Zuerst extrahieren wir die Gennamen.

In der Regex Name=(?P<gene_name>.+?); , +? wird anstelle von + verwendet, weil wir wollen, dass es nicht gierig ist und die Suche beim ersten Semikolon stoppt; Andernfalls stimmt das Ergebnis bis zum letzten Semikolon überein.

Außerdem wird die Regex zuerst mit re.compile kompiliert, anstatt direkt wie in re.search verwendet zu werden, um eine bessere Leistung zu erzielen, da wir sie auf Tausende von Attributzeichenfolgen anwenden werden.

extract_gene_name dient als Hilfsfunktion, die in pd.apply verwendet werden kann. Dies ist die zu verwendende Methode, wenn eine Funktion auf jeden Eintrag eines Datenrahmens oder einer Reihe angewendet werden muss.

In diesem speziellen Fall möchten wir den Gennamen für jeden Eintrag in ndf.attributes und die Namen in einer neuen Spalte namens gene_name wieder zu ndf hinzufügen.

Gen-IDs und Beschreibung werden auf ähnliche Weise extrahiert.

 RE_GENE_ID = re.compile(r'gene_id=(?P<gene_id>ENSG.+?);') 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<desc>.+?);') 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)

Die Regex für RE_GENE_ID ist etwas spezifischer, da wir wissen, dass jede gene_id mit ENSG beginnen muss, wobei ENS für ensembl und G für Gen steht.

Für Einträge ohne Beschreibung geben wir einen leeren String zurück. Nachdem alles extrahiert wurde, verwenden wir die Attribute-Spalte nicht mehr, also lassen wir sie fallen, um die Dinge mit der Methode .drop schön und sauber zu halten:

 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...

Im obigen Aufruf gibt attributes die spezifische Spalte an, die wir löschen möchten.

axis=1 bedeutet, dass wir eine Spalte anstelle einer Zeile löschen ( standardmäßig axis=0 ).

inplace=True bedeutet, dass das Löschen auf dem DataFrame selbst ausgeführt wird, anstatt eine neue Kopie mit der angegebenen gelöschten Spalte zurückzugeben.

Ein kurzer .head Blick zeigt, dass die Attribute-Spalte tatsächlich weg ist und drei neue Spalten hinzugefügt wurden: gene_name , gene_id und desc .

Lassen Sie uns aus Neugier sehen, ob alle gene_id und gene_name eindeutig sind:

 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,)

Überraschenderweise ist die Anzahl der Gennamen kleiner als die der Gen-IDs, was darauf hinweist, dass einige gene_name mehreren Gen-IDs entsprechen müssen. Lassen Sie uns herausfinden, was sie sind.

 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

Wir gruppieren alle Einträge nach dem Wert von gene_name und zählen dann die Anzahl der Elemente in jeder Gruppe mit .count() .

Wenn Sie die Ausgabe von ndf.groupby('gene_name').count() , werden alle Spalten für jede Gruppe gezählt, aber die meisten haben dieselben Werte.

Beachten Sie, dass NA-Werte beim Zählen nicht berücksichtigt werden, nehmen Sie also nur die Anzahl der ersten Spalte, seqid (wir verwenden .ix[:, 0] , um sicherzustellen, dass es keine NA-Werte gibt).

Sortieren Sie dann die Zählwerte mit .sort_values ​​und kehren Sie die Reihenfolge mit .ix[::-1] um.

Als Ergebnis kann ein Genname mit bis zu sieben Gen-IDs geteilt werden.

 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]

Ein genauerer Blick auf alle SCARNA20-Gene zeigt, dass sie tatsächlich alle unterschiedlich sind.

Obwohl sie den gleichen Namen haben, befinden sie sich an unterschiedlichen Positionen des Genoms.

Ihre Beschreibungen scheinen jedoch nicht sehr hilfreich zu sein, um sie zu unterscheiden.

Der Punkt hier ist zu wissen, dass Gennamen nicht für alle Gen-IDs eindeutig sind und etwa 0,15 % der Namen von mehreren Genen geteilt werden.

Wie lang ist ein typisches Gen?

Ähnlich wie wir es getan haben, als wir versuchten, die Unvollständigkeit des Genoms zu verstehen, können wir einfach eine length zu ndf hinzufügen:

 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() berechnet einige einfache Statistiken basierend auf den Längenwerten:

• Die mittlere Länge eines Gens beträgt etwa 36.000 Basen

• Die mittlere Länge eines Gens beträgt etwa 5.200 Basen

• Die minimale und maximale Genlänge beträgt etwa acht bzw. 2,3 Millionen Basen.

Da der Mittelwert viel größer als der Median ist, impliziert dies, dass die Längenverteilung rechtsschief ist. Um einen konkreteren Blick darauf zu werfen, lassen Sie uns die Verteilung grafisch darstellen.

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

Pandas bietet eine einfache Schnittstelle zu matplotlib, um das Plotten mit DataFrames oder Serien sehr praktisch zu machen.

In diesem Fall heißt es, dass wir ein Histogramm ( kind='hist' ) mit 50 Bins wollen und die y-Achse auf einer logarithmischen Skala liegen lassen ( logy=True ).

Aus dem Histogramm können wir sehen, dass sich die Mehrheit der Gene im ersten Bin befindet. Einige Genlängen können jedoch mehr als zwei Millionen Basen betragen. Lassen Sie uns herausfinden, was sie sind:

 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

Wie Sie sehen können, heißt das längste Gen CNTNAP2, was die Abkürzung für Contactin Associated Protein-like 2 ist. Laut seiner Wikipedia-Seite

Dieses Gen umfasst fast 1,6 % des Chromosoms 7 und ist eines der größten Gene im menschlichen Genom.

In der Tat! Wir haben das gerade selbst verifiziert. Was ist dagegen mit den kleinsten Genen? Es stellt sich heraus, dass sie bis zu acht Basen kurz sein können.

 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

Die Längen der beiden Extremfälle liegen fünf Größenordnungen auseinander (2,3 Millionen vs. 8), was enorm ist und ein Hinweis auf die Vielfalt des Lebens sein kann.

Ein einzelnes Gen kann über einen Prozess namens alternatives Spleißen in viele verschiedene Proteine ​​übersetzt werden, etwas, das wir noch nicht erforscht haben. Solche Informationen befinden sich auch in der GFF3-Datei, aber außerhalb des Rahmens dieses Beitrags.

Genverteilung zwischen Chromosomen

Als letztes möchte ich auf die Genverteilung zwischen den Chromosomen eingehen, die auch als Beispiel für die Einführung der .merge Methode zum Kombinieren zweier DataFrames dient. Intuitiv beherbergen längere Chromosomen wahrscheinlich mehr Gene. Mal sehen, ob das stimmt.

 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

Wir haben die chrs Variable aus dem vorherigen Abschnitt ausgeliehen und sie verwendet, um die unassemblierten Sequenzen herauszufiltern. Gemessen an der Ausgabe hat das größte Chromosom 1 tatsächlich die meisten Gene. Während Chromosom Y die kleinste Anzahl von Genen hat, ist es nicht das kleinste Chromosom.

Beachten Sie, dass es im Mitochondrium (MT) anscheinend keine Gene gibt, was nicht stimmt.

Eine etwas stärkere Filterung des ersten von pd.read_csv zurückgegebenen pd.read_csv - df zeigt, dass alle MT-Gene aus der Quelle insdc stammen (die zuvor beim Generieren von edf herausgefiltert wurden, wo wir nur Quellen von havana, ensembl oder ensembl_havana berücksichtigten).

 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...

This example also shows how to combine two conditions during filtering with & ; the logical operator for “or” would be | .

Note that the parentheses around each condition are required, and this part of the syntax in Pandas is different from Python, which would have been literal and and or .

Next, let's borrow the gdf DataFrame from the previous section as a source for the length of each chromosome:

 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

The columns that are not relevant to the analysis are dropped for clarity.

Yes, .drop can also take a list of columns and drop them altogether in one operation.

Note that the row with a seqid of MT is still there; we will get back to it later. The next operation we will perform is merge the two datasets based on the values of 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

Since chr_gene_counts is still a Series object, which doesn't support a merge operation, we need to convert it to a DataFrame object first with .to_frame .

.reset_index() converts the original index (ie seqid ) into a new column and resets current index as 0-based incremental numbers.

The output from cdf.head(2) shows what it looks like. Next, we used the .merge method to combine the two DataFrame on the seqid column ( on='seqid' ).

After merging gdf and cdf , the MT entry is still missing. This is because, by default, .merge operates an inner join, while left join, right join, or outer join are available by tuning the how parameter.

Please refer to the documentation for more details.

Later, you may find that there is also a related .join method. .merge and .join are similar yet have different APIs.

According to the official documentation says

The related DataFrame.join method, uses merge internally for the index-on-index and index-on-column(s) joins, but joins on indexes by default rather than trying to join on common columns (the default behavior for merge). If you are joining on index, you may wish to use DataFrame.join to save yourself some typing.

Basically, .merge is more general-purpose and used by .join .

Finally, we are ready to calculate the correlation between chromosome length and 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

By default .corr calculates the Pearson correlation between all pairs of columns in a dataframe.

But we have only a single pair of columns in this case, and the correlation turns out to be positive – 0.73.

In other words, the larger the chromosome, the more likely it is to have more genes. Let's also plot the two columns after sorting the value pairs by length :

 ax = merged[['length', 'gene_count']].sort_values('length').plot(x='length', y='gene_count',) # 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)) 

As seen in image above, even though it is a positive correlation overall, it does not hold for all chromosomes. In particular, for Chromosome 17, 16, 15, 14, 13, the correlation is actually negative, meaning the number of genes on the chromosome decreases as the chromosome size increases.

Findings and Future Research

That ends our tutorial on the manipulation of an annotation file for human genome in GFF3 format with the SciPy stack. The tools we've mainly used include IPython, Pandas, and matplotlib. During the tutorial, not only have we learned some of the most common and useful operations in Pandas, we also answered some very interesting questions about our genome. In summary:

1. Ungefähr 0,37 % des menschlichen Genoms sind immer noch unvollständig, obwohl der erste Entwurf vor über 15 Jahren herauskam.
2. Basierend auf dieser speziellen GFF3-Datei, die wir verwendet haben, gibt es etwa 42.000 Gene im menschlichen Genom.
3. Die Länge eines Gens kann von einigen Dutzend bis zu über zwei Millionen Basen reichen.
4. Gene sind nicht gleichmäßig auf die Chromosomen verteilt. Insgesamt gilt: Je größer das Chromosom, desto mehr Gene beherbergt es, aber für eine Untergruppe der Chromosomen kann die Korrelation negativ sein.

Die GFF3-Datei ist sehr reich an Anmerkungsinformationen, und wir haben gerade erst an der Oberfläche gekratzt. Wenn Sie an weiteren Erkundungen interessiert sind, sind hier ein paar Fragen, mit denen Sie spielen können:

1. Wie viele Transkripte hat ein Gen normalerweise? Wie viel Prozent der Gene haben mehr als 1 Transkript?
2. Wie viele Isoformen hat ein Transkript typischerweise?
3. Wie viele Exons, CDS und UTRs hat ein Transkript normalerweise? Welche Größen haben sie?
4. Ist es möglich, die Gene anhand ihrer Funktion zu kategorisieren, wie in der Beschreibungsspalte beschrieben?