PAC 2

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)

##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

SECCIÓ 1. Fonaments de programació i bases de dades

EXERCICI 1

Una cadena d’ADN és una cadena de caràcters formada per les lletres “A”, “C”, “G” i “T”. Donada una cadena d’ADN corresponent a una cadena codificant, la seva cadena d’ARN transcrita es forma substituint-hi totes les aparicions de “T” a la cadena d’ADN per “U”.

a) Creeu una funció anomenada DNA_to_RNA que faci la conversió d’ADN a ARN transcrit. L’entrada ha de ser una variable de tipus caràcter com a “GATTACA” i la sortida ha de ser una variable de tipus caràcter com a “GAUUACA” (1 punt).

variable_entrada <- "GATTACA"
DNA_to_RNA<-gsub("T","U",variable_entrada)
variable_salida <- DNA_to_RNA
variable_salida

b) Proveu la funció per a les cadenes d’ADN “GATGTAACTTGACTACGTAAATT” i “ACGGTTCCAATGTAA” (0.25 punts).

variable_entrada <- "GATGTAACTTGACTACGTAAATT"

DNA_to_RNA <-function(variable_entrada) {
 funcion<- gsub("T","U", variable_entrada)
  return(funcion)
}
variable_salida <- DNA_to_RNA(variable_entrada1)
variable_salida


variable_entrada <- "ACGGTTCCAATGTAA"

DNA_to_RNA <-function(variable_entrada) {
 funcion<- gsub("T","U", variable_entrada)
  return(funcion)
}
variable_salida <- DNA_to_RNA(variable_entrada)
variable_salida

EXERCICI 2

Donades dues cadenes d’ADN, es pot definir una noció de distància comptant el nombre de símbols que difereixen entre elles. Per exemple, si tenim les cadenes ACGGT i ACGAC, la distància és 2.

a) Creeu una funció anomenada DNA_distance que calculi la distància entre dues cadenes d’ADN. L’entrada han de ser dues cadenes d’ADN (de la mateixa longitud) i la sortida ha de ser la distància (1 punt).

install.packages("devtools")
devtools::install_github("hiraethus/Needleman-Wunsch")

generar_sequencia_adn <- function(longitud) {
  bases <- c("A","T","C","G")
  sequencia <- sample(bases, longitud,     replace=TRUE)
  return(paste(sequencia, collapse = "")) 
}

seq1 <- generar_sequencia_adn (8)
seq2 <- generar_sequencia_adn (8)
seq1
seq2

library("hiraethus.needleman.wunsch")

match <- 1
mismatch <- -1
gap <- -1

alignment <- needle(seq1, seq2, match)
print(alignment)

b) Proveu la funció per als parells de cadenes d’ADN següents (0.25 punts): “GATGTAACTTGACTACGTAAATT” i “GATGTAACTTGACTGCGTAAATT”; “ACGGTTCCAATGTAA” i “AGGATTCCAATGTAG”

seq3 <- "GATGTAACTTGACTACGTAAATT"
seq4 <- "GATGTAACTTGACTGCGTAAATT"

match <- 1
mismatch <- -1
gap <- -1

alignment2 <- needle(seq3, seq4, match, mismatch, gap)
print(alignment2)

seq5 <- "ACGGTTCCAATGTAA"
seq6 <- "AGGATTCCAATGTAG"

match <- 1
mismatch <- -1
gap <- -1

alignment3 <- needle(seq5, seq6, match, mismatch, gap)
print(alignment3)

EXERCICI 3

En aquest exercici, treballarem amb el conjunt de dades gapminder, que a R es troba al paquet library(gapminder).

Podeu instal·lar el paquet amb la instrucció install.packages(“gapminder”). Un cop tingueu el paquet instal·lat, podreu carregar les dades amb la instrucció data(gapminder).

install.packages("gapminder")
library(gapminder)
data(gapminder)
View(gapminder)

Les dades contenen les variables següents: ● País: País ● Continent: Continent on és el país ● Year: Any de la recollida de les dades ● lifeExp: Esperança de vida ● pop: Població del país ● gdpPercap: PIB per càpita, en dòlars nord-americans, ajustats per inflació

Resoleu les preguntes següents:

a) Amb la funció sqldf, creeu una taula amb els següents requeriments (0.75 punts):

i) Seleccioneu només dades de 2007

ii) Seleccioneu només els països el PIB per càpita dels quals el 2007 sigui menor que 800$ o major que 35000$

iii) Ordeneu les dades de manera ascendent pel PIB per càpita

library(sqldf)
head(sqldf("SELECT * FROM gapminder where year =2007 AND gdpPercap<800 OR gdpPercap>35000 ORDER BY gdpPercap ASC"))

b) Comenteu els resultats de la taula anterior. En particular, us demanem que expliqueu el que veieu sobre els continents, l’esperança de vida i el PIB per càpita. Adjunteu un gràfic que mostri la relació entre PIB per càpita i l’esperança de vida i interpreteu-lo (0.5 punts).

Al indicar en la funció que mostre sols els països que en l’any 2007 tenien una renda per càpita menor que 800 o major que 35000 de manera ascendent, apareixen com a resultat tan sols els països que en aquest any indicat tenien una renda menro que 800 ya que es la instrucció indica primer. Els països amb la menor renda per càpita son tots països d’Àfrica, a més correspon que a un major valor del PIB major son los anys d’esperança de vida, independement de la població. El PIB per càpita i l’esperança de vida presenten una correlació directa.

plot(gapminder$lifeExp, gapminder$gdpPercap, type='h', xlab="Esperanza de vida", ylab="PIB")

En aquest gràfic observem com en països on l’esperança de vida es més elevada, son els que tenen una renda per càpita més elevada encara que hi ha algunes excepcions on l’esperança de vida no es molt elevada, i la renda per capita es exageradament elevada.

c) Amb la funció dbGetQuery, creeu una taula amb els següents requeriments (0.75 punts):

i) Seleccioneu només dades de 2007

ii) Calculeu la mitjana d’esperança de vida i PIB per càpita per continent`

install.packages("DBI")
library(DBI)
library(RSQLite)
library(dplyr)
library(ggplot2)


db <- (gapminder)
db


dades2 <- dbGetQuery("SELECT * FROM gapminder, AVG(lifeExp) AS mitjana_vida, AVG(gdpPercap) AS mitjana_PIB WHERE year = 2007")
dades2

d) Comenteu els resultats de la taula anterior. Torneu a fer un gràfic que relacioni l’esperança de vida i el PIB per càpita. Relacioneu el que veieu aquí amb la vostra anàlisi a la part b) (0.5 punts).

SECCIÓ 2. Probabilitat i simulació

EXERCICI 4

La distribució de Weibull és útil a l’anàlisi de supervivència, que és una àrea de la bioestadística que es dedica a estudiar, per exemple, el temps fins que un pacient es recupera o mor d’una malaltia. A R, la distribució de Weibull està implementada en les funcions dweibull (densitat), pweibull (probabilitat acumulada), qweibull (quantils) i rweibull (simulació de dades aleatòries). Us recomanem que llegiu la documentació de R d’aquestes distribucions (la podeu veure amb la instrucció ?dweibull).

Suposeu que el temps de vida en anys d’un cert tipus de rates de laboratori, en condicions normals, es pot estimar bé per una distribució de Weibull amb paràmetre de forma (shape) és igual a 5 i el paràmetre d’escala (scale) igual a 3. Després d’interferir amb l’expressió d’un gen, la distribució del temps de vida en anys es pot aproximar bé per una distribució de Weibull amb paràmetre igual a 6 i paràmetre d’escala igual a 2.

a) Grafiqueu les densitats dels temps de vida per als ratolins en condicions normals i per als ratolins als quals se’ls modifica l’expressió d’un gen. Comenteu els resultats (0.6 punts).

curve(dweibull(x, shape=5, scale=3), col="blue")
curve(dweibull(x, shape=6, scale=2), col="red", add=TRUE)

Al modificar el gen, el temps de vida en anys de les rates es veu també afectat respecte a el temps de vida de una rata en condicions normals. Una rata a la qual se li ha modificat el gen, viurà menys temps que una a la que no se li ha modificat, ja que com es veu en el gràfic la corba roja, que indica el temps de vida d’aquelles rates modificades, desapareix abans que la corba blava de les rates normals.

b) Calculeu les probabilitats que el temps de vida sigui més gran que dos anys per als dos tipus de ratolins. Després, calculeu la probabilitat que el temps de vida sigui inferior a un any per a ambdós tipus de ratolins. Feu els càlculs de manera exacta i comenteu els resultats (0.6 punts).

pweibull(2, shape=5, scale=3, lower.tail=FALSE, log.p=FALSE)
pweibull(2, shape=6, scale=2, lower.tail=FALSE, log.p=FALSE)
pweibull(1, shape=5, scale=3, lower.tail=TRUE, log.p=FALSE)
pweibull(1, shape=6, scale=2, lower.tail=TRUE, log.p=FALSE)

En aquest cas, els resultats segueixen un poc la dinàmica del apartat anterior, en que es mostra que quan els ratolins son modificats es més probable que visquen menys d’1 any, i menys probable que visquen més de 2 anys. El contrari es reflexa en el cas dels ratolins normals, ja que aquest solen tindre més temps de vida que els modificats genèticament.

c) A través de simulacions, estimeu la mitjana, la variància i els percentils 25%, 50% i 75% del temps de vida per a ambdós tipus de ratolins. Compareu els resultats dels percentils obtinguts a través de simulacions amb els càlculs exactes obtinguts amb la funció qweibull. Feu 10000 simulacions per a cada tipus de ratolí (0.8 punts).

mean(sample(5:3))
var(sample(5:3))

mean(sample(6:2))
var(sample(6:2))

qweibull(0.25,shape=5, scale=3)
qweibull(0.50,shape=5, scale=3)
qweibull(0.75,shape=5, scale=3)
qweibull(0.25,shape=6, scale=2)
qweibull(0.50,shape=6, scale=2)
qweibull(0.75,shape=6, scale=2)

d) Redacteu una breu reflexió sobre els aspectes ètics de fer servir animals en experiments científics (0.5 punts).

Gràcies als comités ètics l’ús d’animals d’experimentació cada vegada es més restringida.La funció principal d’aquests órgans es reemplaçar, reduir i refinar l’ús dels animals en els experiments científics amb l’objectiu de que els científics siguen més conscients sobre el que comporten i en fajen un ús responsable. Aquesta tasca del comité ha estat questionada ja que a l’hora d’iniciar un projecte d’investigació, aquest ha de passar abans pels comites ètics i revisar-ho, però el que ocorre es que de vegades els animals que permeten utilizar son pocs per a que els resultats de l’estudi siguen significatius. Tot i això pense que aquests comites fan una gran tasca per permitir l’ús d’aquests animals i es veu com realment volen obtindre benefici del seu ús però de manera responsable sense arribar a causar-los un gran dany. També crec que encara ens queda molt per avançar en ciència, i encara es veu com una cosa impossible, el no treballar amb animals, ja que moltes vegades es l’objecte d’estudi del cual s’han d’extraure certes coses per a poder estudiar, al igual que no és el mateix utilizar certes línies celulares inmortalizades que fer l’estudi sobre persones. En conclusió, pense que hem d’utilitzar animals d’experimentació en molts projectes sempre que no hi haja una altra alternativa, i si s’han d’utilitzar sempre fer-ho de manera conscient, per produir el menor mal possible a l’animal.

EXERCICI 5

En aquest exercici, calcularem probabilitats relacionades amb una prova diagnòstica per a una malaltia amb càlculs exactes i simulacions. La malaltia que ens interessa afecta l’1% de la població. Si un individu té la malaltia, la probabilitat que la prova diagnòstica sigui positiva és del 98% (i la probabilitat que sigui negativa és del 2%). Si un individu no té la malaltia, la probabilitat que la prova diagnòstica sigui positiva és del 5% (i la probabilitat que sigui negativa és del 95%). La distribució binomial us pot ser útil per resoldre l’exercici.

a) Si escollim 100 individus a l’atzar, quina és la probabilitat que hi hagi almenys un individu que tingui la malaltia? Calculeu aquesta probabilitat mitjançant un càlcul exacte i un càlcul basat en 10000 simulacions. Compareu els resultats (0.6 punts).

n=100 #nombre de individus
p=0.01 #probabilitat que el individu tinga la malaltia
val1=1 #valor exacte
dbinom(val1,n,p)

n=10000
dbinom(val1,n,p)

Es pot observar com al aumentar el nombre d’observacions, disminueix la probabilitat de poder patir la malaltia.

b) Si tenim 10 individus malalts, quina és la probabilitat que totes les proves diagnòstiques siguin positives? Si tenim 990 individus sans, quina és la probabilitat que totes les proves diagnòstiques siguin negatives? Calculeu les probabilitats mitjançant càlculs exactes i 10000 simulacions. Compareu els resultats (0.6 punts).

La sensibilitat és la probabilitat que la prova diagnòstica sigui positiva condicionat al fet que la persona estigui malalta.

L’especificitat és la probabilitat condicionada que la prova diagnòstica sigui negativa atès que la persona no està malalta.

El que interessa saber es la probabilitat condicionada que quan la prova diagnòstica siga positiva la persona estigui realment malalta, és el que es coneix com Valor Predictiu Positiu (VPP). Mentre que la probabilitat que quan la prova diagnòstica sigui negativa la persona estigui realment sana, és el concepte de Valor Predictiu Negatiu (VPN).

sensibilitat <- 0.98
especificitat <- 0.95
prevalença <- 0.01

vpp <- (sensibilitat*prevalença)/(sensibilitat*prevalença+(1-especificitat)*(1-prevalença))
vpp

n=10
p=0.1652614
p10 <- dbinom (10,n,p)
p10

sensibilitat <- 0.98
especificitat <- 0.95
prevalença <- 0.01

vpn <- (especificitat*(1-prevalença))/(especificitat*(1-prevalença)+(1-sensibilitat)*prevalença)
vpn

n=900
p=0.9997874
p900 <- dbinom (900,n,p)
p900

La probabilitat de que 10 persones malaltes tinguen un resultat positiu en les proves es de 1.519547e-08, el que es una probabilitat molt baixeta de que just en 10 persones malaltes totes les proves siguen positives. En canvi, per als individus sans, la probabilitat de que la prova diagnòstica siga negativa es una probabilitat elevada, d’un 0.825835. Això no significa que la prova siga més o menys vàlida, ja que per saber això ens deuriem fixar unicament en el valor predictiu positiu o negatiu corresponent. El que ocorre en aquest cas es que els individus malalts es una n menuda i això fa que el valor predictiu positiu siga un valor molt menut, mentre que els individus sans es una n molt elevada i el valor predictiu negatiu aumenta la fiabilitat del test.

c) Suposem que tenim 10 individus malalts i 990 individus sans. Grafiqueu la distribució de probabilitat del nombre de proves diagnòstiques positives per als 10 individus malalts i els 990 individus sans. Comenteu els resultats (0.7 punts).

plot(dbinom(0:1000,10,0.98))

plot(dbinom(0:1000,900,0.5))

d) Basat en els resultats dels apartats anteriors i el fet que la malaltia afecta a l’1% de la població, si observem una prova diagnòstica positiva, justifiqueu quina opció us sembla més probable (0.6 punts):

i) És més probable que la persona estigui malalta que no ho estigui.

ii) És més probable que la persona no estigui malalta.

iii) No ho sabem.

La opció més probable es la ii) És més probable que la persona no estigui malalta. Ja que com hem vist el valor predictiu positiu de la prova diagnòstica no és del tot fiable, ja que ix una probabilitat molt baixa. De manera que si ix un valor predictiu positiu, el més segur es que la persona no estigui malalta, encara que hi ha una xicoteta probabilitat de que pot ser si ho estiga.