Εισαγωγή στον προγραμματισμό με τη γλώσσα python, Project εξαμήνου

python_project.ipynb

Επαναδιατύπωση του "Αφαιρέστε όσα SNPs έχουν Hardy Weinberg Equilibrium < 0.001"

H συγκεκριμένη έφραση δεν είναι σωστή. Το σωστό είναι: "Αφαιρέστε όσα SNPs έχουν Hardy Weinberg expectation p-value < 0.001"

H wikipedia δίνει ένα πολύ καλό παράδειγμα:
https://en.wikipedia.org/wiki/Hardy%E2%80%93Weinberg_principle#Significance_tests_for_deviation

Πρoσέξτε πως υπολογίζει τα expected values (1467.4, 141.2, 3.4) με βάση τα observed values (1469, 138, 5)

Αφού τα υπολογίσετε μετά μπορείτε να εφαρμόσετε το chi square test στη python:

>>> from scipy import stats
>>> s = stats.chisquare([1469, 138, 5], [1467.4, 141.2, 3.4])
>>> print (s)
Power_divergenceResult(statistic=0.82720700518383861, pvalue=0.66126308839474002)

Προσέξτε ότι το statistic (0.8272..) είναι το ίδιο που δίνει και η wikipedia (0.83), απλά είναι στρογγυλοποιημένο.
Οπότε αυτό που πρέπει να πάρετε είναι το pvalue:

>>> print (s.pvalue)
0.66126308839474002

Αν αυτό το pvalue είναι μικρότερο του 0.001 τότε "πετάμε" το SNP

ΠΡΟΣΟΧΗ: Αυτός είναι ένας από τους πολλούς τρόπους για να το υπολογίσετε! Συγκεκριμένα είναι το asymptotic test (δείτε και: http://pngu.mgh.harvard.edu/~purcell/plink/summary.shtml#hardy). Διάφοροι άλλοι τρόποι μπορεί να δώσουν (λίγο) διαφορετικά αποτελέσματα (φυσικά όλοι οι σωστά υλοποιημένοι τρόποι θα θεωρηθούν σωστοί)

Σημείωση για το "Βρείτε τα statistical significant SNPs (association test)"

Το threshold για significance είναι πολύ μικρότερο από 0.001 !!

Εδώ: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3154648/
στο κεφάλαιο "Multiple testing" έχει εκτενή συζήτηση για το ποιο είναι το σωστό threshold.
Αν θέλετε μπορείτε να εφαρμόσετε "κατά γράμμα" τις οδηγίες που λέει εφαρμόζοντας είτε
Bonferroni είτε Sidak correction.

Αλλά μπορούμε να θεωρήσουμε το 10^-6 ένα safe threshold για τη περίπτωση μας.

Σημείωση για το "Ποια variants είναι σε Linkage Disequilibrium με τα statistical significant SNPs που βρήκατε;"

Περιμένουμε ότι όσο απομακρυνόμαστε από το significance SNP τόσο να μικραίνει και το LD
Οπότε απλά πρέπει να ελέγξετε για LD μόνο τα SNPs τα οποία να είναι "κοντά" στο significant SNP.
Πόσο κοντά όμως;;
Εκεί πρέπει να πειραματιστείτε. Μπορούμε να εφαρμόσουμε το εξής heuristic:
Αν έχουμε αποκρυνθεί (upstream/downstream) τόσο ώστε στα τελευταία 1000 SNPs δεν έχουμε βρει ούτε ένα SNP σε LD με το significant, τότε σταματάμε την αναζήτηση.

(Technically με αυτόν τον τρόπο ψάχνουμε για cis LD)

Feel free να εφαρμόσετε κάποιο δικό σας heuristic!

Παράδειγμα Manhattan plot με ένα χρωμόσωμα

Υπάρχει εδώ: http://2.bp.blogspot.com/-w1h3x3rpOLY/TbXVEHCr1cI/AAAAAAAALqI/wiQqpr1HJ6c/s400/2011-04-19+GGD+manhattan4.png
(Αγνοήστε τα χρώματα αν θέλετε)

ΠΡΟΣΟΧΗ!! Λάθος στον υπολογισμό του LD

Στο link που σας έχω δώσεις για τον υπολογισμό του LD: http://rannala.org/books/CUPChap3.pdf
Έχει ένα λάθος!
Στη σελίδα 68.

το λάθος είναι:

Kudos στον Παύλο Παυλίδη που το εντόπισε!!

Για την ολοκλήρωση του project στέλνω τη λίστα με τα SNPs τα οποία χρησιμοποιήθηκαν για την δημιουργία του artificial dataset:

rs910873 32635433 1 1.7 3 
rs6010620 61780283 1 1.2 1.4
rs1015362 32202273 1 3.0 3.0
rs4911414 32193105 0 3.0 3.0
rs2273535 54394948 1 1.5 1.8
rs1064039 23566427 0 1.5 3.0

Η δεύτερη στήλη είναι το locus του SNP. H τρίτη στήλη (1 ή 0) είναι αν το causative allele είναι το reference (0) ή το alternative (1). H τρίτη και τέταρτη στήλη είναι το heterozygous risk και το homozygous αντίστοιχα. Όπως περιγράφεται και εδώ: https://mathgen.stats.ox.ac.uk/genetics_software/hapgen/hapgen2.html

For example, -dl 1085679 1 1.5 2.25 2190692 0 2 4 specifies two disease SNPs, at positions 1085679 and 2190692, and with heterozgyote risks 1.5 and 2, homozygote risks 2.25 and 4, and risk alleles set to 1 and 0 at each SNP respectively

Σε αυτό το gist: https://gist.github.com/kantale/87b4e1b24abf9be85e44fcd7f8324639 υπάρχει η δική μου (Αλέξανδρος) υλοποίηση καθώς και ένα jupyter notebook με οδηγίες για το πως μπορεί να χρησιμοποιηθεί η υλοποίηση για όλα τα tasks του project.

Σαν τελική παρατήρηση να σημειώσω ότι έχω κάνει ένα μεθοδολογικό λάθος στη δημιουργία του artificial dataset. Συγκεκριμμένα υπάρχουν 3 SNPs (rs910873, rs4911414, rs1015362) τα οποία βρίσκονται κοντά το ένα με το άλλο και έχουν εξαιρετικά μεγάλο effect (3.0). Το hapgen2 όμως θεωρεί ότι τα effects που δημιουργούν τα causative SNPs είναι independent! Οπότε είναι εξαιρετικά σπάνιο να υπάρχουν στη πραγματικότητα 3 SNPs τόσο κοντά το ένα με το άλλο με τόσο μεγάλο effect τα οποία να είναι independent. Αυτό "μπέρδεψε" λίγο το hapgen2 και δημιούργησε πολλά statistical significant SNPs (~78). Επίσης αυτό φάνηκε και στο "περίεργο" qqplot. Συγγνώμη αν μπέρδεψε κάποιους αυτό!