sadII_zadanie_zaliczeniowe

2330 days ago by agorska

load("data/ABAG_1.Rdata") 
       
#zapisywanie PROJEKTU do pliku save.image("data/ABAG_1.Rdata") 
       
#install.packages("fAsianOptions") 
       
library(rtracklayer) library(lattice) library(stats) library(graphics) library(sfsmisc) library(seqinr) library(fAsianOptions) 
       
Loading required package: RCurl
Loading required package: bitops





Loading required package: timeDate
Loading required package: timeSeries
Loading required package: fBasics
Loading required package: MASS

Attaching package: ‘fBasics’

The following object(s) are masked from ‘package:base’:

    norm

Loading required package: fOptions
#sciąganie danych robimy tylko raz session <- browserSession("UCSC") query <- ucscTableQuery(session) tnames<-trackNames(query) 
       
gene.query <- ucscTableQuery(session, "knownGene") #head(sort(trackNames(gene.query))) 
       
tableName(gene.query)<-"knownGene" knownGene <- getTable(gene.query) 
       
#tx - transkrypcyjne #cds - coding #exon #surowe dane head(knownGene) 
       
        name chrom strand txStart txEnd cdsStart cdsEnd exonCount   
exonStarts
1 uc001aaa.3  chr1      +   11873 14409    11873  11873         3   
11873,12612,13220,
2 uc010nxr.1  chr1      +   11873 14409    11873  11873         3   
11873,12645,13220,
3 uc010nxq.1  chr1      +   11873 14409    12189  13639         3   
11873,12594,13402,
4 uc009vis.3  chr1      -   14361 16765    14361  14361         4   
14361,14969,15795,16606,
5 uc009vjc.1  chr1      -   16857 17751    16857  16857         2   
16857,17232,
6 uc009vjd.2  chr1      -   15795 18061    15795  15795         5
15795,16606,16857,17232,17605,
                        exonEnds proteinID    alignID
1             12227,12721,14409,           uc001aaa.3
2             12227,12697,14409,           uc010nxr.1
3             12227,12721,14409,    B7ZGX9 uc010nxq.1
4       14829,15038,15942,16765,           uc009vis.3
5                   17055,17751,           uc009vjc.1
6 15947,16765,17055,17368,18061,           uc009vjd.2
#nie wiem czemu ale mamy 60 chromosomów, trzeba będzie sprawdzić czy są różne czy takie same chromy<-unique(knownGene$chrom) chromy 
       
 [1] chr1                  chr2                  chr3               
chr4                 
 [5] chr5                  chr6                  chr7               
chr8                 
 [9] chr9                  chrM                  chrX               
chrY                 
[13] chr10                 chr11                 chr12              
chr13                
[17] chr14                 chr15                 chr16              
chr17                
[21] chr18                 chr19                 chr20              
chr21                
[25] chr22                 chr6_apd_hap1         chr6_cox_hap2      
chr6_dbb_hap3        
[29] chr6_mcf_hap5         chr6_qbl_hap6         chr4_ctg9_hap1     
chr6_mann_hap4       
[33] chr6_ssto_hap7        chrUn_gl000211        chrUn_gl000212     
chrUn_gl000213       
[37] chrUn_gl000214        chrUn_gl000218        chrUn_gl000219     
chrUn_gl000220       
[41] chrUn_gl000221        chrUn_gl000222        chrUn_gl000223     
chrUn_gl000227       
[45] chrUn_gl000228        chrUn_gl000229        chrUn_gl000237     
chrUn_gl000241       
[49] chrUn_gl000243        chrUn_gl000247        chr17_ctg5_hap1    
chr1_gl000191_random 
[53] chr1_gl000192_random  chr4_gl000193_random 
chr4_gl000194_random  chr7_gl000195_random 
[57] chr9_gl000201_random  chr17_gl000204_random
chr17_gl000205_random chr19_gl000209_random
60 Levels: chr1 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 ...
chrY
# usuwamy wszystkie niepotrzebne chromosomy dim(knownGene) dobreNazwy <- c("chr1", "chr2", "chr3", "chr4", "chr5", "chr6", "chr7", "chr8", "chr9", "chr10", "chr11", "chr12", "chr13", "chr14", "chr15", "chr16", "chr17", "chr18", "chr19", "chr20", "chr21", "chr22", "chrX", "chrY") dobreIndeksy <- knownGene$chrom %in% dobreNazwy knownGeneND <- knownGene[dobreIndeksy,] head(knownGeneND) chromy<-unique(knownGeneND$chrom) chromy 
       
[1] 80922    12



        name chrom strand txStart txEnd cdsStart cdsEnd exonCount   
exonStarts
1 uc001aaa.3  chr1      +   11873 14409    11873  11873         3   
11873,12612,13220,
2 uc010nxr.1  chr1      +   11873 14409    11873  11873         3   
11873,12645,13220,
3 uc010nxq.1  chr1      +   11873 14409    12189  13639         3   
11873,12594,13402,
4 uc009vis.3  chr1      -   14361 16765    14361  14361         4   
14361,14969,15795,16606,
5 uc009vjc.1  chr1      -   16857 17751    16857  16857         2   
16857,17232,
6 uc009vjd.2  chr1      -   15795 18061    15795  15795         5
15795,16606,16857,17232,17605,
                        exonEnds proteinID    alignID
1             12227,12721,14409,           uc001aaa.3
2             12227,12697,14409,           uc010nxr.1
3             12227,12721,14409,    B7ZGX9 uc010nxq.1
4       14829,15038,15942,16765,           uc009vis.3
5                   17055,17751,           uc009vjc.1
6 15947,16765,17055,17368,18061,           uc009vjd.2

 [1] chr1  chr2  chr3  chr4  chr5  chr6  chr7  chr8  chr9  chrX 
chrY  chr10 chr11 chr12 chr13 chr14
[17] chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22
60 Levels: chr1 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 ...
chrY
#musimy oczyscic nasze dane z alternatywnrgo splicingu tmp<-cbind(knownGene$txStart,knownGene$txEnd) knownGeneND<-knownGene[! duplicated(tmp),] 
       
#sortowanie knownGeneNDSort <- knownGeneND[order(knownGeneND$chrom,knownGeneND$txStart),] 
       
#funkcja usuwająca nachodzące na siebie geny usunNachodzaceGeny <- function(dane){ licznik = 1 wynik <-rep(0,length(dane$txStart)) while (licznik<=length(dane$txStart)){ maxEnd = licznik if(licznik!=length(dane$txStart)){ while(dane$txStart[licznik+1] == dane$txStart[licznik]){ if(dane$txEnd[licznik+1]>dane$txEnd[maxEnd]){ maxEnd = licznik+1 } licznik = licznik+1 } } wynik[maxEnd] = 1 licznik = licznik+1 } licznik = 1 while (licznik<=length(dane$txStart)){ minStart = licznik if(licznik!=length(dane$txStart)){ while(dane$txEnd[licznik+1] == dane$txEnd[licznik]){ if(dane$txStart[licznik+1]<dane$txStart[minStart]){ minStart = licznik+1 } licznik = licznik+1 } } wynik[minStart] = wynik[minStart] + 1 licznik = licznik+1 } wynik == 2 } 
       
polaczNachodzace <- function(dane, chromosom){ #wyluskiwanie odpowiedniego chromosomu indeksy<-dane$chrom==chromosom geny<- dane[indeksy,] #lista poczatkow i koncow c odp chromosomu licznik = 1 zaczety <- FALSE tmpStart <- 0 tmpEnd <- 0 wynik <-data.frame(txStart=c(), txEnd=c()) while (licznik<=length(geny$txStart)){ if (zaczety){ if (geny$txStart[licznik] <= tmpEnd) { tmpEnd <- geny$txEnd[licznik] } else { wynik <- rbind(wynik, c(tmpStart, tmpEnd)) zaczety <- FALSE } } else{ zaczety <- TRUE tmpStart <- geny$txStart[licznik] tmpEnd <- geny$txEnd[licznik] } licznik <- licznik + 1 } colnames(wynik)<-c("txStart","txEnd") wynik } 
       
#polacz nachodzace dla wszystkich chromosomow, zloz odpowiedni daraframe wynik<-data.frame() #colnames(wynik)<-c("chrom","txStart","txEnd") for(i in chromy){ tmp<-polaczNachodzace(knownGeneNDSort, i) chr<- c(rep(i, length(tmp$txStart))) tmp<-cbind(chr, tmp) wynik<-rbind(wynik, tmp) } 
       
colnames(wynik)<-c("chrom","txStart","txEnd") GotoweGeny<-wynik head(GotoweGeny) tail(GotoweGeny) 
       
  chrom txStart  txEnd
1  chr1   11873  29370
2  chr1   69090  70008
3  chr1  137840 139281
4  chr1  321145 321207
5  chr1  323891 325896
6  chr1  367658 368597
      chrom  txStart    txEnd
17733 chr22 51017386 51021428
17734 chr22 51039130 51049979
17735 chr22 51065018 51066601
17736 chr22 51135950 51160865
17737 chr22 51176651 51183727
17738 chr22 51195513 51222087
#funkcja tworząca posortowaną listę odległości między genami listaDlugosci<-function(dane,chromosom){ indeksy<-dane$chrom==chromosom dane<- dane[indeksy,] wynik<-data.frame() for(i in 1:(length(dane$txStart)-1)){ wynik<-rbind(wynik,max(0, (dane$txStart[i+1] -dane$txEnd[i])) ) } #posortowany ik<-wynik[order(wynik),] #zwracam wektor: is.vector(ik) } 
       
dane<-listaDlugosci(GotoweGeny,"chr1") head(dane) 
       
[1]  2  3  3  7  7 18
head(lepszeGeny) tail(lepszeGeny) 
       
Error in head(lepszeGeny) : 
  error in evaluating the argument 'x' in selecting a method for
function 'head': Error: object 'lepszeGeny' not found
Error in tail(lepszeGeny) : 
  error in evaluating the argument 'x' in selecting a method for
function 'tail': Error: object 'lepszeGeny' not found
#nie zgubiliśmy żadnych genów przy tej naszej operacji! Dobra ilość! unique(GotoweGeny$chrom) 
       
 [1] chr1  chr2  chr3  chr4  chr5  chr6  chr7  chr8  chr9  chrX 
chrY  chr10 chr11 chr12 chr13 chr14
[17] chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22
24 Levels: chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chrX chrY
chr10 chr11 chr12 chr13 ... chr22
#funkcja zliczająca odległości genów od siebie działa też dla nowego zestawu ! daneDlugosci<-listaDlugosci(GotoweGeny,"chr1") 
       
# head od długosci między genam dla GotoweGeny wygląd dobrze head(daneDlugosci,20) 
       
 [1]  2  3  3  7  7 18 24 24 36 47 55 59 60 69 72 73 73 74 79 83
#czy działa dla pojedynczego ? temp<-listaDlugosci(GotoweGeny,"chr2") length(temp) m<- rep(NA, max_len - length(temp)) temp<-rbind(temp,m) length(temp) class(temp) 
       
[1] 1145
Error: object 'max_len' not found
Error in rbind(temp, m) : object 'm' not found
[1] 1145
[1] "numeric"
#mamy listę list odległości -> kolejne elementy nazywane jak odleglosci<-list() for (i in chromy){ odleglosci[i] <- list(listaDlugosci(GotoweGeny,i)) } 
       
summary(odleglosci) #okej mamy listę odleglosci dla kazdego chromosomu 
       
      Length Class  Mode   
chr1  1731   -none- numeric
chr2  1145   -none- numeric
chr3   918   -none- numeric
chr4   665   -none- numeric
chr5   734   -none- numeric
chr6   984   -none- numeric
chr7   881   -none- numeric
chr8   588   -none- numeric
chr9   733   -none- numeric
chrX   700   -none- numeric
chrY    85   -none- numeric
chr10  684   -none- numeric
chr11 1049   -none- numeric
chr12  832   -none- numeric
chr13  310   -none- numeric
chr14  622   -none- numeric
chr15  925   -none- numeric
chr16  703   -none- numeric
chr17  993   -none- numeric
chr18  266   -none- numeric
chr19 1110   -none- numeric
chr20  432   -none- numeric
chr21  211   -none- numeric
chr22  413   -none- numeric
head(odleglosci$chr1,10) tail(odleglosci$chr1,10) 
       
 [1]  2  3  3  7  7 18 24 24 36 47
 [1]  1464076  1515333  1521071  1530798  1743296  1843462  2692959 
3044220  3108683 21306958
par(mfrow=c(3,2)) licznik=1 for (i in odleglosci){ hist( i , breaks=100, main = paste("Histogram of ", chromy[licznik], sep="")) licznik = licznik+1 } dev.off() 
       
null device 
          1 



par(mfrow=c(3,2)) licznik =1 for (i in odleglosci){ hist(log(i), breaks=100, main = paste("Histogram of" , chromy[licznik]), xlab="log(liczba genów dla przedziału odległości)", ylab="log(częstotliwość)") licznik = licznik+1 } dev.off() 
       
null device 
          1 



#dla większości chromosomów rozkład jest mniej więcej podobny, chociaż dla chromosomu 15 jest bardziej wyrównany. #chromosom Y ma inny rozkład -> to logiczne i sensowne bo chromosom Y jest dużo mniejszy od pozostałych chromosomów i niesie tylko kilkanaście genów 
       
#trzeba zamienić na matrix żeby zrobić box plot a żeby to zrobić trzeba uzupełnić NA do maksymalnej długości 
       
m_odleglosci[3] summary(odleglosci) 
       
Error: object 'm_odleglosci' not found
      Length Class  Mode   
chr1  1731   -none- numeric
chr2  1145   -none- numeric
chr3   918   -none- numeric
chr4   665   -none- numeric
chr5   734   -none- numeric
chr6   984   -none- numeric
chr7   881   -none- numeric
chr8   588   -none- numeric
chr9   733   -none- numeric
chrX   700   -none- numeric
chrY    85   -none- numeric
chr10  684   -none- numeric
chr11 1049   -none- numeric
chr12  832   -none- numeric
chr13  310   -none- numeric
chr14  622   -none- numeric
chr15  925   -none- numeric
chr16  703   -none- numeric
chr17  993   -none- numeric
chr18  266   -none- numeric
chr19 1110   -none- numeric
chr20  432   -none- numeric
chr21  211   -none- numeric
chr22  413   -none- numeric
#tak się nie da bo różne długości wektorów as.matrix(odleglosci$ch1) 
       
Error in array(x, c(length(x), 1L), if (!is.null(names(x)))
list(names(x),  : 
  attempt to set an attribute on NULL
#BoxPloty -> będziemy mogli porównać rozkłady z różnych z chromosomów boxplot(lapply(odleglosci, log), horizontal = TRUE, col=c("blue","red"), notch=TRUE) dev.off() 
       
null device 
          1 
# wygląda na to że rozkład odległości jest w pewnym zakresie podobny. 
       
#QQploty par(mfrow=c(3,2)) licznik =1 for (i in odleglosci){ qqnorm(i, xlab = "Teoretyczne", ylab = "Moje dane", main = paste("QQplot of" , chromy[licznik])) licznik = licznik+1 } dev.off() 
       
null device 
          1 



#QQploty par(mfrow=c(3,2)) licznik =1 for (i in odleglosci){ qqnorm(log(i), xlab = "Teoretyczne", ylab = "Moje dane", main = paste("QQplot of" , chromy[licznik])) licznik = licznik+1 } dev.off() # log(dane) są normalne! 
       
null device 
          1 



sumaDlugosci <- function(chromosom) { suma <- 0 indeksy<-GotoweGeny$chrom==chromosom dane<- GotoweGeny[indeksy,] for(i in 1:(length(dane$txStart)-1)){ suma <- suma + max(0, (dane$txStart[i+1] - dane$txEnd[i])) } as.numeric(suma) } class(sumaDlugosci("chrY")) # suma długości działa dobrze 
       
[1] "numeric"
#robimy data frame z sumami odległości między genami SumaOdleglosci <- data.frame(chr=c(),dlugosci=c()) for (i in chromy){ tmp<-data.frame(chr=i, dlugosci=as.numeric(sumaDlugosci(i))) SumaOdleglosci <-rbind(SumaOdleglosci,tmp ) } colnames(SumaOdleglosci) <- c("chr", "dlugosci") 
       
class(SumaOdleglosci$dlugosci) summary(SumaOdleglosci$dlugosci) 
       
[1] "numeric"
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
 22660000  56180000  95920000  91910000 123400000 177300000 
plot(x=SumaOdleglosci$chr, y=as.vector(SumaOdleglosci$dlugosci),"b", main= "Sumy odleglosci miedzy genami", sub="dla poszczegolnych chromosomow", col=c("blue","red")) dev.off() # ten wykres wyglądał kiedyś zupełnie inaczej 
       
null device 
          1 
# chcę sprawdzić jak ro się ma do wielkości chromosomu: session <- browserSession("UCSC") query <- ucscTableQuery(session, "cytoBand") #tableNames(query) tableName(query)<-"cytoBand" cband <- getTable(query) 
       
chr.sizes<-with(cband, tapply(chromEnd, chrom,max)) barchart(~sort(chr.sizes),col=c("blue","red")) dev.off() 
       
null device 
          1 
#tworzę data.frame z nazwą chromosomu, wielkością chromosomu i sumą odległości między genami v.odleglosci<-as.vector(SumaOdleglosci$dlugosci) suma.a.wielkosc <-cbind(chr.sizes, v.odleglosci) suma.a.wielkosc <- as.data.frame(suma.a.wielkosc) colnames(suma.a.wielkosc) <- c("chr.size", "dlugosci") suma.a.wielkosc class(suma.a.wielkosc) 
       
       chr.size  dlugosci
chr1  249250621 177261491
chr10 135534747 175233100
chr11 135006516 135458025
chr12 133851895 145762212
chr13 115169878 134750497
chr14 107349540 122418171
chr15 102531392 112279224
chr16  90354753 110289551
chr17  81195210 109966562
chr18  78077248 126327462
chr19  59128983  56369010
chr2  243199373  97408776
chr20  63025520  97355032
chr21  48129895  94486234
chr22  51304566  77041866
chr3  198022430  61695086
chr4  191154276  55627050
chr5  180915260  66533074
chr6  171115067  54463382
chr7  159138663  58941141
chr8  146364022  38969918
chr9  141213431  44746173
chrX  155270560  29881262
chrY   59373566  22656136
[1] "data.frame"
#nazwy chromosomów w kolejności zgodnej z kolejnością w data.frame suma.a.wielkosc chr.chrsize <- row.names(suma.a.wielkosc) 
       
xyplot(dlugosci~chr.size, data=suma.a.wielkosc) dev.off() 
       
null device 
          1 
cor(suma.a.wielkosc$chr.size, suma.a.wielkosc$dlugosci) 
       
[1] 0.05759699
#hm.. wygląda jakby była jakas prawidłowość #teraz spróbuję dopasować odległości do długości chromosomów #albo może tu jest błąd? Albo tu -> robi się z błędem fit<-lm(dlugosci~chr.size, data=suma.a.wielkosc ) fit 
       
Call:
lm(formula = dlugosci ~ chr.size, data = suma.a.wielkosc)

Coefficients:
(Intercept)     chr.size  
  8.626e+07    4.381e-02  
#analiza dopasowania modelu liniowego, wykresy mówią że nawet udało się to dopasować! layout(matrix(1:4,2,2)) plot(fit) dev.off() 
       
null device 
          1 
%html Tam mamy wewnątrz tego tracka "Gene bounds", co chyba bedzie najlepsze do zlokalizowania gdzie są geny. Mamy jeszcze "Assembly", tam będzie cała sekwencja wszystkiego i z tego weźmiemy CG. 
       
Tam mamy wewnątrz tego tracka "Gene bounds", co chyba bedzie najlepsze do zlokalizowania gdzie są geny. Mamy jeszcze "Assembly", tam będzie cała sekwencja wszystkiego i z tego weźmiemy CG.
# parametr gęstości genów na kawałku o koordynatach poczatek i koniec lambdaG<- function(chromosom, poczatek, koniec){ lambdaG <- 0 ostatnia <- poczatek indeksy <- GotoweGeny$chrom==chromosom dane <- GotoweGeny[indeksy,] #ustalić który gen zaczyna się jest na poczatku licznik <- 1 while(TRUE){ if(dane$txStart[licznik]>=poczatek){ break }else{licznik = licznik+1}} while((ostatnia < (koniec-1)) && (licznik <= length(dane$txStart))){ ostatnia = dane$txStart[licznik] lambdaG = lambdaG+1 licznik = licznik+1 } lambdaG } lambdaG("chr2" , 1909000, 7000000) 
       
[1] 14
#prostszy zliczacz genów: lambdaG.caly.chrom<- function(chromosom){ length(odleglosci$chromosom) } lambdaG.caly.chrom("chr1") 
       
[1] 0
indeksy<-GotoweGeny$chrom=="chr2" dane<- GotoweGeny[indeksy,] tail(dane) 
       
     chrom   txStart     txEnd
2873  chr2 242615158 242617529
2874  chr2 242641455 242668896
2875  chr2 242674029 242708231
2876  chr2 242749406 242758739
2877  chr2 242792032 242801058
2878  chr2 242912833 242919427
chr.chrsize[2] #GotoweGeny[1] indeksy<-GotoweGeny$chrom==chr.chrsize[2] dane<- GotoweGeny[indeksy,] lambda=length(t(dane)) lambda 
       
[1] "chr10"




[1] 2055
#lista lambdG -> gęstości genów na chromosomach lambdaG.lista <- data.frame(chr=c(), lambda=c()) licznik = 1 while (licznik <= length(chr.chrsize)){ indeksy<-GotoweGeny$chrom==chr.chrsize[licznik] dane<- GotoweGeny[indeksy,] lambda <- length(t(dane)) tmp<-data.frame(chr=chr.chrsize[licznik],lambda) lambdaG.lista <-rbind(lambdaG.lista, tmp) licznik = licznik+1 } suma.a.wielkosc <- cbind(suma.a.wielkosc, lambdaG.lista) 
       
lambdaG.lista 
       
     chr lambda
1   chr1   5196
2  chr10   2055
3  chr11   3150
4  chr12   2499
5  chr13    933
6  chr14   1869
7  chr15   2778
8  chr16   2112
9  chr17   2982
10 chr18    801
11 chr19   3333
12  chr2   3438
13 chr20   1299
14 chr21    636
15 chr22   1242
16  chr3   2757
17  chr4   1998
18  chr5   2205
19  chr6   2955
20  chr7   2646
21  chr8   1767
22  chr9   2202
23  chrX   2103
24  chrY    258
suma.a.wielkosc 
       
       chr.size  dlugosci   chr lambda
chr1  249250621 177261491  chr1   5196
chr10 135534747 175233100 chr10   2055
chr11 135006516 135458025 chr11   3150
chr12 133851895 145762212 chr12   2499
chr13 115169878 134750497 chr13    933
chr14 107349540 122418171 chr14   1869
chr15 102531392 112279224 chr15   2778
chr16  90354753 110289551 chr16   2112
chr17  81195210 109966562 chr17   2982
chr18  78077248 126327462 chr18    801
chr19  59128983  56369010 chr19   3333
chr2  243199373  97408776  chr2   3438
chr20  63025520  97355032 chr20   1299
chr21  48129895  94486234 chr21    636
chr22  51304566  77041866 chr22   1242
chr3  198022430  61695086  chr3   2757
chr4  191154276  55627050  chr4   1998
chr5  180915260  66533074  chr5   2205
chr6  171115067  54463382  chr6   2955
chr7  159138663  58941141  chr7   2646
chr8  146364022  38969918  chr8   1767
chr9  141213431  44746173  chr9   2202
chrX  155270560  29881262  chrX   2103
chrY   59373566  22656136  chrY    258
#For these reasons, we have developed a method that accounts for variation in the distance between consecutive genes, gene density and GC content to identify extreme chromosomal regions of gene expression. 
       
#p,r koordynaty na chromosomie, ale sumaDlugosciKoordynaty<-function(chromosom,pocz,koniec){ indeksy<-GotoweGeny$chrom==chromosom dane<- GotoweGeny[indeksy,] suma <- 0 #ustalić który gen zaczyna się jest na poczatku licznik <- 1 while(TRUE) { if(dane$txStart[licznik]>=pocz){ break } else{ licznik = licznik+1 } } poprzednia <- pocz while(poprzednia < (koniec-1)){ poprzednia = dane$txEnd[licznik] if(dane$txStart[licznik+1]<koniec){ suma <- suma + max(0, (dane$txStart[licznik+1] - dane$txEnd[licznik])) } licznik = licznik+1 } as.numeric(suma) } sumaDlugosciKoordynaty("chr1",11873,249153315) 
       
[1] 177246360
integrand <- function(u){(102*u)^(1)*exp((-1)*102*u)} w <- integrate(integrand, lower = 0, upper = 520998) w class(w) 
       
0 with absolute error < 0
[1] "integrate"
# Specifically, the probability of observing a cluster of r + 1 genes over a base pair distance as short or shorter than the observed value of Si,r is computed from the density f(si,r) as follows: #tutaj nie działa całka # r - ilość interwałów pomiędzy genami #prawdopodobieństwo zaobserwowania klastra na obszarze r+1 genów ja obszarze krótszym niż Sir to : prawdopodobienstwoKlastra<- function(GotoweGeny, chromosom, i,r){ indeksy<-GotoweGeny$chrom==chromosom dane<- GotoweGeny[indeksy,] poczatek <- dane$txStart[i] koniec <- dane$txEnd[i+r] lambda <- lambdaG(chromosom, poczatek, koniec) gammaR <- gamma(r) S <- sumaDlugosciKoordynaty(chromosom, poczatek, koniec) #calka integrand <- function(u){(lambda*u)^(r-1)*exp((-1)*lambda*u)} w <- integrate(integrand, lower = 0, upper = S) if(w) { wynik <- lambda/gammaR * as.numeric(w)} else{ return 0 } wynik } prawdopodobienstwoKlastra(GotoweGeny, "chr1", 1,10) 
       
WARNING: Output truncated!  
full_output.txt





















Error: unexpected numeric constant in:
"        else{
        return 0"
Error: unexpected '}' in "        }"
      chrom   txStart     txEnd
1      chr1     11873     29370
2      chr1     69090     70008
3      chr1    137840    139281
4      chr1    321145    321207
5      chr1    323891    325896
6      chr1    367658    368597
7      chr1    566092    566115
8      chr1    566239    566263
9      chr1    621095    622034
10     chr1    661138    668479
11     chr1    671823    671885
12     chr1    674239    679736
13     chr1    761585    762902
14     chr1    763063    789726
15     chr1    846814    850328
16     chr1    860529    894679
17     chr1    896828    899229
18     chr1    910578    917473
19     chr1    948846    949919
20     chr1    995082   1001833
21     chr1   1007196   1007946
22     chr1   1017197   1051736
23     chr1   1102483   1102578
24     chr1   1104384   1104467
25     chr1   1109285   1121243
26     chr1   1138887   1142089
27     chr1   1146719   1149548
28     chr1   1152287   1159348
29     chr1   1177825   1182102
30     chr1   1189291   1209234
31     chr1   1217488   1227409
32     chr1   1227763   1260046
33     chr1   1266725   1269844
34     chr1   1270657   1277504
35     chr1   1288070   1297157
36     chr1   1309109   1310818

...

17683 chr22  45725524  45737836
17684 chr22  45765811  45802746
17685 chr22  45898718  45998697
17686 chr22  46020410  46020512
17687 chr22  46067916  46156478
17688 chr22  46316247  46373008
17689 chr22  46435786  46440748
17690 chr22  46449725  46454040
17691 chr22  46481876  46487006
17692 chr22  46508631  46508653
17693 chr22  46509565  46509648
17694 chr22  46546730  46639653
17695 chr22  46651559  46659219
17696 chr22  46691039  46692557
17697 chr22  46711928  46726707
17698 chr22  46726771  46730190
17699 chr22  46731297  46753237
17700 chr22  46764755  46805171
17701 chr22  47033736  47075688
17702 chr22  47080306  47134152
17703 chr22  47159638  47311916
17704 chr22  47857047  47882860
17705 chr22  48020216  48027318
17706 chr22  48027451  48251349
17707 chr22  48885287  48943157
17708 chr22  49176106  49176165
17709 chr22  49808173  50051190
17710 chr22  50166936  50173958
17711 chr22  50277310  50280826
17712 chr22  50312282  50321186
17713 chr22  50356294  50357720
17714 chr22  50497819  50524358
17715 chr22  50528684  50580618
17716 chr22  50587817  50597686
17717 chr22  50624358  50638027
17718 chr22  50644745  50656045
17719 chr22  50656117  50664618
17720 chr22  50683612  50686901
17721 chr22  50688491  50700089
17722 chr22  50703223  50709340
17723 chr22  50713407  50746001
17724 chr22  50781745  50862865
17725 chr22  50883430  50913464
17726 chr22  50920152  50924866
17727 chr22  50941375  50946135
17728 chr22  50946644  50962840
17729 chr22  50964181  50968514
17730 chr22  50968837  50971008
17731 chr22  50989540  51001328
17732 chr22  51007289  51010882
17733 chr22  51017386  51021428
17734 chr22  51039130  51049979
17735 chr22  51065018  51066601
17736 chr22  51135950  51160865
17737 chr22  51176651  51183727
17738 chr22  51195513  51222087
Error: unexpected '}' in "}"
[1] "integrate"
Error in prawdopodobienstwoKlastra(GotoweGeny, "chr1", 1, 10) : 
  (list) object cannot be coerced to type 'double'
#The GC percent track shows the percentage of G (guanine) and C (cytosine) bases in 5-base windows. High GC content is typically associated with gene-rich areas. This track may be configured in a variety of ways to highlight different aspects of the displayed information. Click the "Graph configuration help" 
       
#pobieranie informacji o profilach GC library(IRanges) session <- browserSession("UCSC") query <- ucscTableQuery(session) tnames<-trackNames(query) gene.query <- ucscTableQuery(session, "gc5Base") tableName(gene.query) GCpercent <- ucscTableQuery(session, "gc5Base", GRangesForUCSCGenome("hg19", "chr12",IRanges(57795963, 57815592))) GCper<-track(GCpercent) GCper.df<-as.data.frame(GCper) 
       
Attaching package: ‘IRanges’

The following object(s) are masked from ‘package:base’:

    cbind, eval, intersect, Map, mapply, order, paste, pmax,
pmax.int, pmin, pmin.int,
    rbind, rep.int, setdiff, table, union





NULL
head(GCper.df) 
       
  space    start      end width   score
1 chr12 57795966 57795970     5 39.3701
2 chr12 57795971 57795975     5 59.8425
3 chr12 57795976 57795980     5 39.3701
4 chr12 57795981 57795985     5 79.5276
5 chr12 57795986 57795990     5 39.3701
6 chr12 57795991 57795995     5 59.8425
zawartoscGC<-function(chr, pocz, koniec){ session <- browserSession("UCSC") query <- ucscTableQuery(session) tnames<-trackNames(query) gene.query <- ucscTableQuery(session, "gc5Base") tableName(gene.query) GCpercent <- ucscTableQuery(session, "gc5Base", GRangesForUCSCGenome("hg19", chr, IRanges(pocz, koniec))) GCper<-track(GCpercent) GCper.df<-as.data.frame(GCper) mean(GCper.df$score) } zawartoscGC("chr12", 57795963, 57815592) 
       
[1] 42.16019
# teraz mamy już wszystko