text:
function(d, y.lim=NA, new=T, pcol=c(0,1,1), ns = .01){
#Compute mean and 95% confidence interval for each column of d
#and plot as a ribbon plot

	nitems <- dim(d)[2]
	o <- matrix(0,3,nitems)
	d <- d + matrix(runif(prod(dim(d))), dim(d)[1], dim(d)[2]) * ns
	for(i1 in c(1:nitems)){
		t <- t.test(d[,i1])
		o[1,i1] <- t$estimate
		o[2:3,i1] <- t$conf.int[1:2]
		}
	if(is.na(y.lim[1])) y.lim <- range(o)
	if(new) plot(0,0, xlim = c(1, nitems), ylim=y.lim, ylab="Accuracy", xlab="Time")
	lines(1:nitems, o[1,], type = "o", pch=16, lwd=2, col = hsv(pcol[1], pcol[2], pcol[3]))
	polygon(c(1:nitems, nitems:1), c(o[2,], o[3,nitems:1]), col=hsv(pcol[1],pcol[2],pcol[3],.2), border=NA)
	invisible(o)
}

ext: r


text:
# Goal: The amazing R vector notation.

cat("EXAMPLE 1: sin(x) for a vector --\n")
# Suppose you have a vector x --
x = c(0.1,0.6,1.0,1.5)

# The bad way --
n = length(x)
r = numeric(n)
for (i in 1:n) {
  r[i] = sin(x[i])
}
print(r)

# The good way -- don't use loops --
print(sin(x))

cat("\n\nEXAMPLE 2: Compute the mean of every row of a matrix --\n")
# Here's another example. It isn't really about R; it's about thinking in
# matrix notation. But still.
# Let me setup a matrix --
N=4; M=100;
r = matrix(runif(N*M), N, M)

# So I face a NxM matrix
#               [r11 r12 ... r1N]
#               [r21 r22 ... r2N]
#               [r32 r32 ... r3N]
# My goal: each column needs to be reduced to a mean.

# Method 1 uses loops:
mean1 = numeric(M)
for (i in 1:M) {
  mean1[i] = mean(r[,i])
}

# Alternatively, just say:
mean2 = rep(1/N, N) %*% r               # Pretty!

# The two answers are the same --
all.equal(mean1,mean2[,])
#
# As an aside, I should say that you can do this directly by using
# the rowMeans() function. But the above is more about pedagogy rather
# than showing you how to get rowmeans.

cat("\n\nEXAMPLE 3: Nelson-Siegel yield curve\n")
# Write this asif you're dealing with scalars --
# Nelson Siegel function
nsz <- function(b0, b1, b2, tau, t) {
  tmp = t/tau
  tmp2 = exp(-tmp)
  return(b0 + ((b1+b2)*(1-tmp2)/(tmp)) - (b2*tmp2))
}

timepoints <- c(0.01,1:5)

# The bad way:
z <- numeric(length(timepoints))
for (i in 1:length(timepoints)) {
  z[i] <- nsz(14.084,-3.4107,0.0015,1.8832,timepoints[i])
}
print(z)

# The R way --
print(z <- nsz(14.084,-3.4107,0.0015,1.8832,timepoints))

cat("\n\nEXAMPLE 3: Making the NPV of a bond--\n")
# You know the bad way - sum over all cashflows, NPVing each.
# Now look at the R way.
C = rep(100, 6)
nsz(14.084,-3.4107,0.0015,1.8832,timepoints)        # Print interest rates
C/((1.05)^timepoints)                               # Print cashflows discounted @ 5%
C/((1 + (0.01*nsz(14.084,-3.4107,0.0015,1.8832,timepoints))^timepoints)) # Using NS instead of 5%
# NPV in two different ways --
C %*% (1 + (0.01*nsz(14.084,-3.4107,0.0015,1.8832,timepoints)))^-timepoints
sum(C * (1 + (0.01*nsz(14.084,-3.4107,0.0015,1.8832,timepoints)))^-timepoints)
# You can drop back to a flat yield curve at 5% easily --
sum(C * 1.05^-timepoints)

# Make a function for NPV --
npv <- function(C, timepoints, r) {
  return(sum(C * (1 + (0.01*r))^-timepoints))
}
npv(C, timepoints, 5)

# Bottom line: Here's how you make the NPV of a bond with cashflows C
# at timepoints timepoints when the zero curve is a Nelson-Siegel curve --
npv(C, timepoints, nsz(14.084,-3.4107,0.0015,1.8832,timepoints))
# Wow!

# ---------------------------------------------------------------------------
# Elegant vector notation is amazingly fast (in addition to being beautiful)
N <- 1e5
x <- runif(N, -3,3)
y <- runif(N)

method1 <- function(x,y) {
  tmp <- NULL
  for (i in 1:N) {
    if (x[i] < 0) tmp <- c(tmp, y[i])
  }
  tmp
}

method2 <- function(x,y) {
  y[x < 0]
}

s1 <- system.time(ans1 <- method1(x,y))
s2 <- system.time(ans2 <- method2(x,y))
all.equal(ans1,ans2)
s1/s2           # On my machine it's 2000x faster

max_stars_repo_path: RFrontEndSolution/Tests Repo/Rscripts All (163 files)/b5.r


text:
library(ggvis)

mtcars %>%
  ggvis(x = ~wt) %>%
  layer_histograms(width = 1) %>%
  save_spec("layer/histogram.json")

mtcars %>%
  ggvis(x = ~wt, stroke = ~cyl) %>%
  group_by(cyl) %>%
  layer_freqpolys(width = 1) %>%
  save_spec("layer/freqpoly-grouped.json")

mtcars %>%
  ggvis(x = ~wt, y = ~mpg) %>%
  layer_smooths() %>%
  save_spec("layer/smooth.json")

mtcars %>%
  ggvis(x = ~wt, y = ~mpg) %>%
  group_by(cyl) %>%
  layer_smooths() %>%
  save_spec("layer/smooth-grouped.json")

max_stars_repo_name: romainfrancois/ggvis

max_issues_repo_name: romainfrancois/ggvis

max_forks_repo_name: romainfrancois/ggvis


text:
library(ltm)

dat <- read.csv("result.csv", header=F)
mod  <- rasch(dat, IRT.param=TRUE)
print(mod)
plot(mod)
print(summary(mod))

max_stars_repo_name: crotsu/ItemResponseTheory

max_issues_repo_name: crotsu/ItemResponseTheory

max_forks_repo_name: crotsu/ItemResponseTheory


text:
Sys.setlocale("LC_ALL", "C")

#library(reshape2)
library(ggplot2)

monthly_stats <- function(DF) {

    #print(DF)
    #DF1 <- melt(DF, id.var="Date")

    #print(DF1)
    plt <- ggplot(DF, aes(x = Date, y = F2)) +
            ylab("效率(%)") +
            #ggtitle(iconv("月运行情况","UTF-8","UTF-8")) + #iconv
            geom_bar(stat = "identity", alpha=0.9, fill="#2ECC40", position="dodge") +
            scale_x_continuous(breaks=DF$Date, labels=DF$DateString) +
            theme(axis.text.x = element_text(angle=90, vjust=0.1, size=6, ,  color="#BBBFC2"),

               axis.title.x=element_blank(),
                    axis.text.y = element_text(size=7,  color="#BBBFC2"),
                    axis.title.y = element_text(size=10,  color="#BBBFC2"),

                                panel.grid.major=element_line(color="#444444"),
                    panel.grid.minor=element_line(color="#444444"),
                    panel.background=element_blank(),
                    #panel.border=element_rect(color= "#1F1F1F"),
                    legend.title = element_text(size=7,  color="#BBBFC2"),
                    legend.text = element_text(size=7,  color="#BBBFC2"),
                    legend.background = element_rect(colour = '#1F1F1F', fill = '#1F1F1F'),
                    plot.background=element_rect(fill = "#1F1F1F", color = "#1F1F1F")

            )

    ggsave("output/monthly_stats002.png",width=6, height=3)

}

DF <- read.table("monthly2.csv", sep=",", header=TRUE)

monthly_stats(DF)


text:
#library(ggplot2)
library(data.table)
library(BiocParallel)
library(limma)
library(edgeR)
library(variancePartition)
library(dplyr)
library(Matrix)

pdf(file = "figures.pdf") # The height of the plot in inches

geneExpr = fread('./output/pyscenic_aucell_output.csv') %>% as.data.frame
rownames(geneExpr) = geneExpr$V1
geneExpr = geneExpr[,-1]
geneExpr = as.matrix(geneExpr)

metadata = fread('[PATH]/metadata.csv')
metadata$is_covid = factor(metadata$is_covid)
metadata$subject_id = factor(metadata$subject_id)
metadata$tissue = factor(metadata$tissue)
metadata$pid = factor(metadata$pid)
metadata$celltype = factor(metadata$celltype)
metadata$Biosample = factor(with(metadata, paste(subject_id, pid, tissue)))
metadata$aggr_unit = factor(with(metadata, paste(subject_id, pid, tissue_id,  celltype))) # tissue_id
metadata$ID = factor(with(metadata, paste(subject_id, pid, tissue, celltype)))

aggregate_counts = function( geneExpr, metadata){
  # aggrecate gene counts
  geneCounts = lapply( levels(metadata$aggr_unit), function(id){
      idx = which(metadata$aggr_unit == id)
      colMeans(geneExpr[idx,,drop=FALSE])
    })
  geneCounts = do.call(cbind, geneCounts)
  colnames(geneCounts) = levels(metadata$aggr_unit)
  # simplify metadata
  md2 = lapply( levels(metadata$aggr_unit), function(id){
      unique(metadata[metadata$aggr_unit == id,c("subject_id", "Biosample", "clinical", "is_covid", "pid", "tissue", "celltype")])
  })
  md2 = do.call(rbind, md2)
  list( geneCounts  = geneCounts,
        metadata    = md2)
}
res = aggregate_counts( geneExpr, metadata)
form = ~ (1|tissue) + (1|celltype) + (1|subject_id)
#vobj = voomWithDreamWeights(res$geneCounts, form, res$meta, plot=TRUE, BPPARAM=SnowParam(24, progressbar=TRUE))

## dream
test_per_celltype = function(CT){
  idx = which(info$celltype == CT)
  form = ~ 0 + is_covid:tissue + (1|tissue) + (1|subject_id)
  L1  = getContrast(mat[,idx], form, info[idx,],
    c("is_covid1:tissuePFC",'is_covid0:tissuePFC'))
  L2  = getContrast(mat[,idx], form, info[idx,],
    c("is_covid1:tissuemedulla",'is_covid0:tissuemedulla'))
  L3  = getContrast(mat[,idx], form, info[idx,],
    c("is_covid1:tissuechoroidplexus",'is_covid0:tissuechoroidplexus'))
  L = data.frame(PFC = L1, medulla = L2, choroidpleus = L3)
  fit = dream(mat[,idx], form, info[idx,], L = L, BPPARAM=SnowParam(24, progressbar=TRUE))
  fit
    }

info = res$metadata
mat = res$geneCounts

fitList = lapply( levels(info$celltype), test_per_celltype)
names(fitList) = levels(info$celltype)

resTab = lapply( names(fitList), function(CT){
  fit = fitList[[CT]]
  tabList = lapply( colnames(coef(fit))[1:3], function(tissue){
    tab = topTable(fit, coef=tissue, number=Inf)

    data.frame( Gene      = rownames(tab),
                Tissue    = tissue,
                CellType  = CT,
                logFC     = tab$logFC,
                P.Value   = tab$P.Value,
                z.std     = tab$z.std)
  })
  do.call(rbind, tabList)
})
resTab = data.table(do.call(rbind, resTab))
dev.off()


text:
library(sf)
library(tmap)
library(raster)

dem = raster('../geonomics/geonomics/data/yosemite_demo/CA_DEM.tif')
tmp = raster('../geonomics/geonomics/data/yosemite_demo/yosemite_lyrs/tmp_1980-2010.tif')

#set the bbox values, because this somehow got lost and screwed up when I hackily
#ripped this data out of the annoying-as-hell ArcLayer package that I downloaded it in
#NOTE: values pulled directly from the file's metadata XML file (which I saved in
#this directory)
dem@extent@xmin = -124.482003
dem@extent@xmax = -114.131170
dem@extent@ymin = 32.528832
dem@extent@ymax = 42.009665

#CRS is already lined up, so just clip, then resample
dem_crop = crop(dem, tmp)
dem_crop_resamp = resample(dem_crop, tmp)

#check the resolutions are the same
res(dem_crop_resamp) == res(tmp)

#then just going to use a linear model to backfill the missing values that were clipped
#out of the original DEM layer because they fall in NV; this is fine because the DEM
#is not being used for analysis at all; rather, it's just being used as the terrain
#onto which the results will be draped
dem_vals = dem_crop_resamp[!is.na(dem_crop_resamp)]
tmp_vals = tmp[!is.na(dem_crop_resamp)]
df = data.frame(alt=dem_vals, tmp=tmp_vals)
mod = lm(alt ~ tmp, data=df)
preds = stats::predict.lm(mod, data.frame(tmp=tmp[is.na(dem_crop_resamp)]))
dem_crop_resamp[is.na(dem_crop_resamp)] = preds

#plot to check
tmap_mode('view')
tm_shape(dem_crop_resamp) + tm_raster(alpha = 0.5) + tm_shape(tmp) + tm_raster(alpha = 0.5)

#now write to new file
writeRaster(dem_crop_resamp, 'yosemite_DEM.tif', 'GTiff', overwrite=T)


text:
library(purrr)

# Download as much as we can from Wu's folder on github.
download_from_github = function(file_name,
    github_data_folder = "https://raw.githubusercontent.com/wxwx1993/PM_COVID/master/Data/"
) {
    download.file(
        paste0(github_data_folder, file_name),
        file.path("data", file_name)
    )
}

download_from_github("census_county_interpolated.csv")
download_from_github("county_base_mortality.txt")
download_from_github("county_old_mortality.txt")
download_from_github("county_014_mortality.txt")
download_from_github("county_1544_mortality.txt")
download_from_github("county_4564_mortality.txt")
download_from_github("county_pm25.csv")
download_from_github("state_policy0410.csv")
download_from_github("temp_seasonal_county.csv")

# This is different from the brfss_county_interpolated on github, which has
# BMI's, not obesity rates.
download.file(
    "https://www.countyhealthrankings.org/sites/default/files/media/document/analytic_data2020.csv",
    file.path("data", "analytic_data2020.csv")
)

download.file(
    "https://opendata.arcgis.com/datasets/6ac5e325468c4cb9b905f1728d6fbf0f_0.csv?outSR=%7B%22latestWkid%22%3A3857%2C%22wkid%22%3A102100%7D",
    file.path("data", "hospitals.csv")
)

# Download data from a specific commit to replicate as closely as possible.
download_from_JHU = function(date,
    commit = "968f9e627140a318cc5876b79fb959298da9a348"
) {
    github_folder = paste0(
        "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/",
        commit,
        "/csse_covid_19_data/csse_covid_19_daily_reports/"
    )
    file_name = paste0(format(date, "%m-%d-%Y"), ".csv")
    download.file(
        paste0(github_folder, file_name),
        file.path("data", commit, file_name))
}

download_from_JHU(as.Date("2020-03-30"))
walk(
    seq(as.Date("2020-03-22"), as.Date("2020-09-07"), by = "days"),
    download_from_JHU
)


text:
<rePP:domainCheck names="asd.com,example.com,my-ex3.com" _wrap=domainCheck/>
<rePP:domainInfo name="my-ex4.com" _wrap=domainInfo/>

ext: rsx

max_stars_repo_name: hiqsol/reclient

max_issues_repo_name: hiqsol/reclient

max_forks_repo_name: hiqsol/reclient


text:
# -------------------------------------
# Color palettes for plotting functions
# -------------------------------------

### from smillie et al
### https://github.com/cssmillie/ulcerative_colitis

# rcolorbrewer "set" palette
set.colors = c(brewer.pal(9, 'Set1'), brewer.pal(7, 'Set2'), brewer.pal(12, 'Set3')[c(3,9,8)], 'violetred4')
set.colors[6] = 'khaki2'
set.colors[8] = 'lightskyblue2'
set.colors = rep(set.colors, 10)

# material.io "heat" palette
material.heat = function(n){
    colorRampPalette(
        c(
        "#283593", # indigo 800
        "#3F51B5", # indigo
        "#2196F3", # blue
        "#00BCD4", # cyan
        "#4CAF50", # green
        "#8BC34A", # light green
        "#CDDC39", # lime
        "#FFEB3B", # yellow
        "#FFC107", # amber
        "#FF9800", # orange
        "#FF5722"  # deep orange
        )
    )(n)
}

# NMF library "heatmap" palette
nmf.colors = c("#D73027", "#FC8D59", "#FEE090", "#FFFFBF", "#E0F3F8", "#91BFDB", "#4575B4")

max_stars_repo_name: goldfinchbio/organoid-paper

max_issues_repo_name: goldfinchbio/organoid-paper

max_forks_repo_name: goldfinchbio/organoid-paper


text:
REBOL [
  Title:   "Make Red & Rebol lexer tests"
	Author:  "Peter W A Wood"
	File: 	 %make-lexer-auto-tests.r
	Version: 0.1.0
	Tabs:	 4
	Rights:  "Copyright (C) 2016 Peter W A Wood. All rights reserved."
	License: "BSD-3 - https://github.com/red/red/blob/origin/BSD-3-License.txt"
]

;; initialisations
red-tests: copy ""                          ;; string to hold tests
rebol-tests: copy ""
make-dir %auto-tests/
red-file-out: %auto-tests/lexer-auto-test.red
rebol-file-out: %auto-tests/lexer-auto-test.r

;; create Red test file with header
append red-tests "Red [^(0A)"
append red-tests {  Title:   "Red auto-generated lexer red-tests"^(0A)}
append red-tests {	Author:  "Peter W A Wood"^(0A)}
append red-tests {  File: 	 %lexer-auto-test.red^(0A)}
append red-tests {  License: "BSD-3 - https://github.com/dockimbel/Red/blob/origin/BSD-3-License.txt"^(0A)}
append red-tests "]^(0A)^(0A)"
append red-tests "^(0A)^(0A)comment {"
append red-tests "  This file is generated by make-lexer-auto-red-tests.r^(0A)"
append red-tests "  Do not edit this file directly.^(0A)"
append red-tests "}^(0A)^(0A)"
append red-tests join ";make-length:"
                  [length? read %make-lexer-auto-tests.r "^(0A)^(0A)"]
append red-tests {#include  %../../../../quick-test/quick-test.red^(0A)^(0A)}

;; create Red test file with header
append rebol-tests "Rebol [^(0A)"
append rebol-tests {  Title:   "Red auto-generated lexer rebol-tests"^(0A)}
append rebol-tests {	Author:  "Peter W A Wood"^(0A)}
append rebol-tests {  File: 	 %lexer-auto-test.r^(0A)}
append rebol-tests {  License: "BSD-3 - https://github.com/dockimbel/Red/blob/origin/BSD-3-License.txt"^(0A)}
append rebol-tests "]^(0A)^(0A)"
append rebol-tests "^(0A)^(0A)comment {"
append rebol-tests "  This file is generated by make-lexer-auto-rebol-tests.r^(0A)"
append rebol-tests "  Do not edit this file directly.^(0A)"
append rebol-tests "}^(0A)^(0A)"
append rebol-tests {do %../../../../lexer.r^(0A)}
append rebol-tests {do %../../../../quick-test/quick-unit-test.r^(0A)^(0A)}

;; add the actual tests and write out the files
the-tests: read %lexer.tests
append red-tests the-tests
replace/all the-tests "load/all" "lexer/process"
replace/all the-tests "--assert [" "--assert [[] "
replace/all the-tests {lexer/process "} {lexer/process "Red []}
replace/all the-tests "lexer/process {" "lexer/process {Red []"
append rebol-tests the-tests
remove/part find rebol-tests ";red" find/tail rebol-tests ";der"
write red-file-out red-tests
write rebol-file-out rebol-tests
print what-dir

print "lexer auto-test files generated"


text:
#define TWO

#include "conf.h"
#include "SysTypes.r"
#include "Types.r"
#include "BalloonTypes.r"
#include "AERegistry.r"
#include "AEUserTermTypes.r"

type 'STH#' as 'STR#';

#include "Common.r"

#include "Two.r"
#include "FiltDefs.r"
#include "audit.r"
#include "credits.r"

#ifdef DEMO
include "RegDemo.rsrc" not 'ckid';
#else
include "RegTwo.rsrc" not 'ckid';
#endif

max_stars_repo_name: officialrafsan/EUDORA

max_issues_repo_name: evilneuro/eudora-mac

max_forks_repo_name: evilneuro/eudora-mac


text:
\docType{data}
\name{dsö}
\alias{dsö}
\title{Dünya Sağlık Örgütü tüberküloz verileri}
\format{7240 satır ve 60 sütundan oluşan bir tibble
\describe{
\item{ülke}{Ülke adı}
\item{iso2}{İki haneli ISO ülke kodu}
\item{iso3}{Üç haneli ISO ülke kodu}
\item{yıl}{Yıl}
\item{yeni_byp_e014}{Yeni vakalar, balgam yayması pozitif (byp), erkek (e), 0 - 14 yaş aralığı (014)}
\item{yeni_byp_e1524}{Yeni vakalar, balgam yayması pozitif (byp), erkek (e), 15 - 24 yaş aralığı (1524)}
\item{yeni_byp_e2534}{Yeni vakalar, balgam yayması pozitif (byp), erkek (e), 25 - 34 yaş aralığı (2534)}
\item{yeni_byp_e3544}{Yeni vakalar, balgam yayması pozitif (byp), erkek (e), 35 - 44 yaş aralığı (3544)}
\item{yeni_byp_e4554}{Yeni vakalar, balgam yayması pozitif (byp), erkek (e), 45 - 54 yaş aralığı (4554)}
\item{yeni_byp_e5564}{Yeni vakalar, balgam yayması pozitif (byp), erkek (e), 55 - 64 yaş aralığı (5564)}
\item{yeni_byp_e65}{Yeni vakalar, balgam yayması pozitif (byp), erkek (e), 65 yaş üstü (65)}
\item{yeni_byp_k014}{Yeni vakalar, balgam yayması pozitif (byp), kadın (k), 0 - 14 yaş aralığı (014)}
\item{yeni_byp_k1524}{Yeni vakalar, balgam yayması pozitif (byp), kadın (k), 15 - 24 yaş aralığı (1524)}
\item{yeni_byp_k2534}{Yeni vakalar, balgam yayması pozitif (byp), kadın (k), 25 - 34 yaş aralığı (2534)}
\item{yeni_byp_k3544}{Yeni vakalar, balgam yayması pozitif (byp), kadın (k), 35 - 44 yaş aralığı (3544)}
\item{yeni_byp_k4554}{Yeni vakalar, balgam yayması pozitif (byp), kadın (k), 45 - 54 yaş aralığı (4554)}
\item{yeni_byp_k5564}{Yeni vakalar, balgam yayması pozitif (byp), kadın (k), 55 - 64 yaş aralığı (5564)}
\item{yeni_byp_k65}{Yeni vakalar, balgam yayması pozitif (byp), kadın (k), 65 yaş üstü (65)}
\item{yeni_byn_e014}{Yeni vakalar, balgam yayması negatif (byn), erkek (e), 0 - 14 yaş aralığı (014)}
\item{yeni_byn_e1524}{Yeni vakalar, balgam yayması negatif (byn), erkek (e), 15 - 24 yaş aralığı (1524)}
\item{yeni_byn_e2534}{Yeni vakalar, balgam yayması negatif (byn), erkek (e), 25 - 34 yaş aralığı (2534)}
\item{yeni_byn_e3544}{Yeni vakalar, balgam yayması negatif (byn), erkek (e), 35 - 44 yaş aralığı (3544)}
\item{yeni_byn_e4554}{Yeni vakalar, balgam yayması negatif (byn), erkek (e), 45 - 54 yaş aralığı (4554)}
\item{yeni_byn_e5564}{Yeni vakalar, balgam yayması negatif (byn), erkek (e), 55 - 64 yaş aralığı (5564)}
\item{yeni_byn_e65}{Yeni vakalar, balgam yayması negatif (byn), erkek (e), 65 yaş üstü (65)}
\item{yeni_byn_k014}{Yeni vakalar, balgam yayması negatif (byn), kadın (k), 0 - 14 yaş aralığı (014)}
\item{yeni_byn_k1524}{Yeni vakalar, balgam yayması negatif (byn), kadın (k), 15 - 24 yaş aralığı (1524)}
\item{yeni_byn_k2534}{Yeni vakalar, balgam yayması negatif (byn), kadın (k), 25 - 34 yaş aralığı (2534)}
\item{yeni_byn_k3544}{Yeni vakalar, balgam yayması negatif (byn), kadın (k), 35 - 44 yaş aralığı (3544)}
\item{yeni_byn_k4554}{Yeni vakalar, balgam yayması negatif (byn), kadın (k), 45 - 54 yaş aralığı (4554)}
\item{yeni_byn_k5564}{Yeni vakalar, balgam yayması negatif (byn), kadın (k), 55 - 64 yaş aralığı (5564)}
\item{yeni_byn_k65}{Yeni vakalar, balgam yayması negatif (byn), kadın (k), 65 yaş üstü (65)}
\item{yeni_ad_e014}{Yeni vakalar, akciğer dışı (ad), erkek (e), 0 - 14 yaş aralığı (014)}
\item{yeni_ad_e1524}{Yeni vakalar, akciğer dışı (ad), erkek (e), 15 - 24 yaş aralığı (1524)}
\item{yeni_ad_e2534}{Yeni vakalar, akciğer dışı (ad), erkek (e), 25 - 34 yaş aralığı (2534)}
\item{yeni_ad_e3544}{Yeni vakalar, akciğer dışı (ad), erkek (e), 35 - 44 yaş aralığı (3544)}
\item{yeni_ad_e4554}{Yeni vakalar, akciğer dışı (ad), erkek (e), 45 - 54 yaş aralığı (4554)}
\item{yeni_ad_e5564}{Yeni vakalar, akciğer dışı (ad), erkek (e), 55 - 64 yaş aralığı (5564)}
\item{yeni_ad_e65}{Yeni vakalar, akciğer dışı (ad), erkek (e), 65 yaş üstü (65)}
\item{yeni_ad_k014}{Yeni vakalar, akciğer dışı (ad), kadın (k), 0 - 14 yaş aralığı (014)}
\item{yeni_ad_k1524}{Yeni vakalar, akciğer dışı (ad), kadın (k), 15 - 24 yaş aralığı (1524)}
\item{yeni_ad_k2534}{Yeni vakalar, akciğer dışı (ad), kadın (k), 25 - 34 yaş aralığı (2534)}
\item{yeni_ad_k3544}{Yeni vakalar, akciğer dışı (ad), kadın (k), 35 - 44 yaş aralığı (3544)}
\item{yeni_ad_k4554}{Yeni vakalar, akciğer dışı (ad), kadın (k), 45 - 54 yaş aralığı (4554)}
\item{yeni_ad_k5564}{Yeni vakalar, akciğer dışı (ad), kadın (k), 55 - 64 yaş aralığı (5564)}
\item{yeni_ad_k65}{Yeni vakalar, akciğer dışı (ad), kadın (k), 65 yaş üstü (65)}
\item{yeni_nüks_e014}{Yeni nükseden vakalar, erkek (e), 0 - 14 yaş aralığı (014)}
\item{yeni_nüks_e1524}{Yeni nükseden vakalar, erkek (e), 15 - 24 yaş aralığı (1524)}
\item{yeni_nüks_e2534}{Yeni nükseden vakalar, erkek (e), 25 - 34 yaş aralığı (2534)}
\item{yeni_nüks_e3544}{Yeni nükseden vakalar, erkek (e), 35 - 44 yaş aralığı (3544)}
\item{yeni_nüks_e4554}{Yeni nükseden vakalar, erkek (e), 45 - 54 yaş aralığı (4554)}
\item{yeni_nüks_e5564}{Yeni nükseden vakalar, erkek (e), 55 - 64 yaş aralığı (5564)}
\item{yeni_nüks_e65}{Yeni nükseden vakalar, erkek (e), 65 yaş üstü (65)}
\item{yeni_nüks_k014}{Yeni nükseden vakalar, kadın (k), 0 - 14 yaş aralığı (014)}
\item{yeni_nüks_k1524}{Yeni nükseden vakalar, kadın (k), 15 - 24 yaş aralığı (1524)}
\item{yeni_nüks_k2534}{Yeni nükseden vakalar, kadın (k), 25 - 34 yaş aralığı (2534)}
\item{yeni_nüks_k3544}{Yeni nükseden vakalar, kadın (k), 35 - 44 yaş aralığı (3544)}
\item{yeni_nüks_k4554}{Yeni nükseden vakalar, kadın (k), 45 - 54 yaş aralığı (4554)}
\item{yeni_nüks_k5564}{Yeni nükseden vakalar, kadın (k), 55 - 64 yaş aralığı (5564)}
\item{yeni_nüks_k65}{Yeni nükseden vakalar, kadın (k), 65 yaş üstü (65)}
}}
\usage{dsö}
\description{Dünya Sağlık Örgütü Küresel Tüberküloz Raporu verilerinin bir alt kümesi.}
\seealso{\code{\link[tidyr]{who}}}
\keyword{datasets}

ext: rd

max_stars_repo_name: botan/veriler

max_issues_repo_name: botan/veriler

max_forks_repo_name: botan/veriler


text: DF2 <- rbind(DF, DF[7:8,, drop=FALSE])

max_stars_repo_name: DeagleGross/Rplugin

max_issues_repo_name: DeagleGross/Rplugin

max_forks_repo_name: DeagleGross/Rplugin


text:
#Function to run the Beer Game for a single time step
#James Paine
#jpaine@mit.edu

PirateBeerGame_Opt = function(AI_Entity = FALSE, AI_Entity_Index=FALSE, AI_parameter_df = FALSE, AI_Order = FALSE,
                          t, Orders,
                          Order_flows, Shipping_flows, OH_Inventory, Backorder,
                          L, L_hat, Production_Request,
                          Parameter_df = FALSE, Integer_Ordering = FALSE, Noisey_Ordering = FALSE,
                          Shipping_Delay = 2, Information_Delay = 2
                          ) {

  if (Parameter_df != FALSE) {

    #Repair blanks or zero values
    Parameter_df[is.na(Parameter_df)] = 0

    #Assign parameter values
    theta = unlist(Parameter_df$theta)
    alpha_s = unlist(Parameter_df$alpha_s)
    beta = unlist(Parameter_df$beta)
    S_prime = unlist(Parameter_df$S_prime)

    #TeamName = as.character(Parameter_sub$Team_Name[1])

  } else {

    theta = rep(0.36,4)
    alpha_s = rep(0.26,4)
    beta = rep(0.34,4)
    S_prime = rep(17,4)

    TeamName = "Default Average Agents"

  }

  if (unname(AI_parameter_df[1] != FALSE)) {

    theta[AI_Entity_Index] = unname(AI_parameter_df["theta"])
    alpha_s[AI_Entity_Index] = unname(AI_parameter_df["alpha_s"])
    beta[AI_Entity_Index] = unname(AI_parameter_df["beta"])
    S_prime[AI_Entity_Index] = unname(AI_parameter_df["S_prime"])

  }

  #####Recieve Inventory and Advance Shipping Delays#####

  #Recieve shipments

  if (t!= 1) {

    OH_Inventory[,t] = OH_Inventory[,t-1] + Shipping_flows[,1,t-1]

    #Advance shipping delays
    Shipping_flows[,1:(Shipping_Delay-1),t] = Shipping_flows[,2:(Shipping_Delay),t-1]

  }

  #####Fill Orders######

  #View Orders
  Incoming_Order = Order_flows[,1,t-1]
  #Add backlog to order pressure (remove to smooth)
  Incoming_Order = Incoming_Order + Backorder[,t-1]

  #Ship what you can
  #PROBLEM: NEED TO ACCOUNT FOR DEMAND FROM BACK ORDER
  Outbound_shipments = pmax(0,pmin(OH_Inventory[,t],Incoming_Order))

  #Put shipments into lefthand shipping slot
  Shipping_flows[1:3,Shipping_Delay,t] = Outbound_shipments[2:4]
  #Shipping_flows[1:3,2,t] = Outbound_shipments[2:4]

  #Send shipments from retailer to the customer
  Final_Customer_Orders_Filled[t] = Outbound_shipments[1]

  #Determine backlog, if any
  Inventory_Shortage = Order_flows[,1,t-1] - OH_Inventory[,t]
  Backorder[,t] = Backorder[,t-1] + Inventory_Shortage
  Backorder[,t][Backorder[,t]<0] = 0

  #Subtract out the order from inventory
  OH_Inventory[,t] = OH_Inventory[,t] - Outbound_shipments

  #####Advance Order Slips and Brewers Brew######

  Order_flows[,1:(Information_Delay-1),t] = Order_flows[,2:Information_Delay,t-1]
  #Order_flows[,1,t] = Order_flows[,2,t-1]
  Order_flows[1,1,t] = Orders[t]

  #Brewers Brew
  Shipping_flows[4,Shipping_Delay,t] = Production_Request[t-1]
  #Shipping_flows[4,2,t] = Production_Request[t-1]

  #####PLACE ORDERS######

  for (i in 1:4) {

    Entity_Index = i

    #Observe supply line
    SL = sum(Shipping_flows[Entity_Index,,t]) #Total supply line inbound

    #L[Entity_Index,t] = Incoming_Order[Entity_Index]
    L[Entity_Index,t] = max(min(Incoming_Order[Entity_Index],OH_Inventory[Entity_Index,t]),Order_flows[Entity_Index,1,t-1])
    L[Entity_Index,t] = Order_flows[Entity_Index,1,t-1]

    if (t==1) {
      L_hat[Entity_Index,t] = Order_flows[Entity_Index,1,t]
    } else if (t==2) {
      L_hat[Entity_Index,t-1] = Order_flows[Entity_Index,1,t-1]
      L[Entity_Index,t-1] = L_hat[Entity_Index,t-1]
    }

    L_hat[Entity_Index,t] = theta[Entity_Index]*L[Entity_Index,t] + (1-theta[Entity_Index])*L_hat[Entity_Index,t-1]

    S = OH_Inventory[Entity_Index,t] #Stock as of current time
	  S = OH_Inventory[Entity_Index,t] - Backorder[Entity_Index,t] #Stock incorporating backorder

    #Add noise to the order if needed
    if (Noisey_Ordering == TRUE) {
      eps = rnorm(1, Noise_Mean, Noise_SD)
    } else {
      eps = 0
    }

    if ((AI_Entity == TRUE) & (Entity_Index == AI_Entity_Index)) { #AI Decision

      ####AI DECISION GOES HERE

      if (AI_Order != FALSE) {

        Order_Placed = AI_Order

      } else {

        Order_Placed = max(0,L_hat[Entity_Index,t]+alpha_s[Entity_Index]*(S_prime[Entity_Index]-S-beta[Entity_Index]*SL) + eps)

      }

      #########################

    } else { #'Human Player' decisions

      Order_Placed = max(0,L_hat[Entity_Index,t]+alpha_s[Entity_Index]*(S_prime[Entity_Index]-S-beta[Entity_Index]*SL) + eps)

    }

    # #TURN ON FOR INTEGER ONLY ORDERING
    if (Integer_Ordering == TRUE) {
      Order_Placed = round(Order_Placed,0)
    }

    if (Entity_Index == 4) {
      Production_Request[t] = Order_Placed
    } else {
      Order_flows[Entity_Index+1,Information_Delay,t] = Order_Placed
    }

  } #next i

  # #Debug
  #
  # Order_flows, Shipping_flows, OH_Inventory, Backorder,
  # L, L_hat, Production_Request
  #
  # str(Order_flows)
  # str(Shipping_flows)
  # str(OH_Inventory)
  # str(Backorder)
  # str(L)
  # str(L_hat)
  # str(Production_Request)

  Final_Orders = unname(rbind(Order_flows[2:4,Information_Delay,],Production_Request))
  Amplification = (abs(Final_Orders[,t]-Orders[t]))/Orders[t]

  Lagged_Orders = vector()
  for (n in 2:4) {
    frontmatter = rep(NA,((n-1)*Information_Delay))
    Lagged_Orders = rbind(Lagged_Orders,append(frontmatter,Orders[-(length(Orders)-(n-1)*Information_Delay+1):-length(Orders)]))

  }
  Lagged_Orders=unname(rbind(Orders,Lagged_Orders))

  LaggedAmplifcation = (abs(Final_Orders[,t]-Lagged_Orders[,t]))/Lagged_Orders[,t]
  TotalLaggedAmpCosts = sum(LaggedAmplifcation^2,na.rm=TRUE)

  #KEY VARIABLE TO MINIMIZE OVER
  Max_Amp = max(Amplification)
  reward = abs(sum(-25*Amplification^2 + 1))
  #reward = TotalLaggedAmpCosts

  fnt_output = list(Order_flows, Shipping_flows, OH_Inventory, Backorder,
                    L, L_hat, Production_Request, Amplification, reward)

  names(fnt_output) = c("Order_flows", "Shipping_flows", "OH_Inventory", "Backorder",
                        "L", "L_hat", "Production_Request","Amplification", "reward")

  return(fnt_output)

} #END FUNCTION


text:
context("tbl_sql")

test_that("can generate sql tbls with raw sql", {
  df <- data_frame(x = 1:3, y = 3:1)
  tbls <- test_load(df, ignore = "df")

  test_f <- function(tbl) {
    clone <- tbl(x$src, build_sql("SELECT * FROM ", x$from))
    expect_equal(collect(tbl), collect(clone))
  }
})

test_that("NAs in character fields handled by db sources (#2256)", {
  df <- data.frame(
    x = c("a", "aa", NA),
    y = c(NA, "b", "bb"),
    z = c("cc", NA, "c"),
    stringsAsFactors = FALSE
  )
  tbls <- test_load(df, ignore = "df")
  for (tbl in tbls) {
    expect_equal(collect(tbl), as.tbl(df))
  }
})


text:
#DPR environment
#ml GCC/6.4.0-2.28 Python/3.6.3 OpenMPI/2.1.1 numpy/1.13.1-Python-3.6.3 pandas/0.18.1-Python-3.6.3 scipy/0.19.1-Python-3.6.3  tabix/0.2.6 R
#source /data/c***/z***/tools/TIGAR/tiger_env/bin/activate
#export PYTHONPATH=$(python -c 'import sys; print(sys.path[-1])'):${PYTHONPATH}

#---R script for DPR model training---
args<-as.numeric(commandArgs(TRUE))

#subjob id
subjob=args

#training tissue
tissue='Brain_Frontal_Cortex_BA9'

#load annotation for gene position
anno<-read.table('/data/c***/z***/anno/gencode/37/gencode.v32.GRCh37.txt',header = T,stringsAsFactors = F)
anno<-anno[,c('chr','left','right','geneid','genename')]

#load expression (residual)
exp_all<-read.table(paste0('/data/c***/z***/data/gtex/exp/v8/exp_residual/',tissue,'.v8.exp.residual'),header = T,stringsAsFactors = F)
colnames(exp_all)<-sapply(colnames(exp_all), function(x) sub('GTEX.','GTEX-',x))
exp<-merge(anno,exp_all,by='geneid')
exp$chr=as.numeric(sub('chr','',exp$chr))
exp<-exp[order(exp$chr),]

#assign 50 genes for each subjob (the 50 genes have to be on the same chromosome)
sub_id=0
for(i in 1:22){
  chr_tmp<-exp[exp$chr==i,1:4]
  chr_times<-ceiling(nrow(chr_tmp)/50)
  chr_tmp$subjob=rep((sub_id+1):(sub_id+chr_times),50)[1:nrow(chr_tmp)]
  sub_id=sub_id+chr_times
  if(i==1){
    chr_df=chr_tmp
  }else{
    chr_df<-rbind(chr_df,chr_tmp)
  }
}
#check DO NOT RUN.
#chr_df<-chr_df[order(chr_df$subjob),]

exp<-exp[which(chr_df$subjob==subjob),]
chr=exp[1,'chr']

#mkdir
dir.create(paste0('/data/g***/z***/gtex/model_training/DPR/raw/',tissue,'/',subjob))
setwd(paste0('/data/g***/z***/gtex/model_training/DPR/raw/',tissue,'/',subjob))
dir.create('./output')

#write expression file (input of DPR)
exp<-exp[,c(2,3,4,1,5,seq(6,ncol(exp)))]
colnames(exp)[1:5]<-c('CHROM','GeneStart','GeneEnd','TargetID','GeneName')
write.table(exp,paste0('./exp.txt'),quote = F,sep='\t',row.names = F)
#write sample file (input of DPR)
sampleid_df<-as.data.frame(colnames(exp)[6:ncol(exp)],stringsAsFactors = F)
write.table(sampleid_df,paste0('./sampleid.txt'),row.names = F,col.names = F,sep='\t',quote = F)

### Train DPR model
data_dir=paste0('/data/g***/z***/gtex/model_training/DPR/raw/',tissue,'/',subjob)

setwd('/data/c***/z***/tools/TIGAR/TIGAR-master/')

cmd=paste0('/data/c***/z***/tools/TIGAR/TIGAR-master/TIGAR_Model_Train.sh --model DPR --Gene_Exp ',data_dir,'/exp.txt --train_sampleID ',data_dir,'/sampleid.txt --genofile /data/c***/z***/data/gtex/geno/v8/liftover/pruning/vcf/gtex_rsid_prune0.9.vcf.gz --chr ',chr,' --genofile_type vcf --Format GT  --maf 0.05 --hwe 0.00001 --cvR2 1 --dpr 1 --ES fixed --out_dir ',data_dir,'/output')
system(cmd,wait = T,ignore.stdout=T,ignore.stderr=T)

max_stars_repo_name: mjbetti/MR-JTI

max_issues_repo_name: mjbetti/MR-JTI

max_forks_repo_name: mjbetti/MR-JTI


text:
context("Basic Tests")

test_that("A basic rgithub context can be acquired", {
  create.github.context("https://api.github.com")
  repos <- get.user.repositories("cscheid")
  repos_overview <- do.call("rbind",
                            lapply(repos$content[1:5], function(x) {
                              data.frame(name = x$name,
                                         owner = x$owner$login,
                                         updated_at = x$updated_at)}))
  cat("\n")
  print(repos_overview)
})

max_stars_repo_name: gordonwoodhull/rgithub

max_issues_repo_name: gordonwoodhull/rgithub

max_forks_repo_name: att/rgithub


text:
#' Retrieve size of content | body | headers
#'
#' @param har_resp_obj HAR response object
#' @family splash_har_helpers
#' @export
get_content_size <- function(har_resp_obj) {
  csize <- har_resp_obj$response$content$size
  if (is.null(csize)) return(NA_real_)
  return(as.numeric(csize))
}

#' @rdname get_content_size
#' @export
get_body_size <- function(har_resp_obj) {
  bsize <- har_resp_obj$response$bodySize
  if (is.null(bsize)) return(NA_real_)
  return(as.numeric(bsize))
}

#' @rdname get_content_size
#' @export
get_headers_size <- function(har_resp_obj) {
  hsize <- har_resp_obj$response$headersSize
  if (is.null(hsize)) return(NA_real_)
  return(as.numeric(hsize))
}

max_stars_repo_name: hrbrmstr/splashr

max_issues_repo_name: hrbrmstr/splashr

max_forks_repo_name: hrbrmstr/splashr


text:
data_path <- "/extraspace/TCGA/TCGA_protein"

file_name <- list.files(data_path, "_protein.re_20160627")
file_name <- file_name[!grepl(file_name, pattern = "GBM|LAML")]

library(tidyverse)
library(stringr)

get_meth <- function(x) {
  print(x)
  name <- str_replace(x, "_.*", "")
  dataset_id <- paste("TCGA", name, sep='-')

  data <- read_tsv(file.path(data_path, x)) %>% select(-length(.))
  type <- str_sub(colnames(data),14,14)
  col_name <- str_sub(colnames(data),1,12)
  col_name[type=="1"] = paste(name,"Normal", col_name[type=="1"], sep="-")
  col_name[type=="0"] = paste(name,"Tumor", col_name[type=="0"], sep="-")
  colnames(data) <- col_name

  data %>%
    subset(select = c(1,which(grepl(name, colnames(.))))) %>%
    subset(select = which(!duplicated(colnames(.)))) %>%
    drop_na()

  snor_pro <- read_tsv(
    file.path(
      "/home/jgong/gongj_snoRNA/database_data/protein",
      list.files(
        path = "/home/jgong/gongj_snoRNA/database_data/protein",
        pattern = paste(name,".*_id", sep = "")
        )
      )
    ) %>%
    mutate(dataset_id = dataset_id) %>%
    select(snorna = sno, dataset_id, protein = gene, spearman_corr = r, p_value=p, fdr)

  snor_pro %>%
    left_join(data, by = "protein") %>%
    gather(sample_id, expression, -c(snorna,dataset_id, protein, spearman_corr, p_value, fdr)) %>%
    drop_na() ->
    result_data
  result_data %>%
    write_tsv(path=paste("/home/cliu18/liucj/projects/4.SNORic/database/protein_expresion/protein_expression", name,"tsv", sep="."),na="", col_names=F)
}

file_name %>%
  lapply(get_meth)

# %>%
  # write_tsv(path = "/home/cliu18/liucj/projects/4.SNORic/database/protein_expresion/protein_expression.tsv", col_names = F)


text:
#create new address field
job_work_merge$address <- paste(job_work_merge$House..,job_work_merge$Street.Name,job_work_merge$Borough,"NY",sep=" ")

# Geocoding a csv column of "addresses" in R

#load ggmap
library(ggmap)

# Select the file from the file chooser
fileToLoad <- file.choose(new = TRUE)

# Read in the CSV data and store it in a variable
origAddress <- read.csv(fileToLoad, stringsAsFactors = FALSE)

# Initialize the data frame
geocoded <- data.frame(stringsAsFactors = FALSE)

# Loop through the addresses to get the latitude and longitude of each address and add it to the
# origAddress data frame in new columns lat and lon
for(i in 1:nrow(origAddress))
{
  # Print("Working...")
  result <- geocode(origAddress$addresses[i], output = "latlona", source = "google")
  origAddress$lon[i] <- as.numeric(result[1])
  origAddress$lat[i] <- as.numeric(result[2])
  origAddress$geoAddress[i] <- as.character(result[3])
}
# Write a CSV file containing origAddress to the working directory
write.csv(origAddress, "geocoded.csv", row.names=FALSE)


text:
\name{bfast-package}
\alias{bfast-package}
\title{Breaks For Additive Season and Trend (BFAST)}

\description{
BFAST integrates the decomposition of time series into trend, seasonal, and remainder
components with methods for detecting and characterizing abrupt changes within the trend and
seasonal components. BFAST can be used to analyze different types of satellite image time
series and can be applied to other disciplines dealing with seasonal or non-seasonal time
series,such as hydrology, climatology, and econometrics. The algorithm can be extended to
label detected changes with information on the parameters of the fitted piecewise linear
models.

Additionally monitoring disturbances in BFAST-type models at the end of time series
(i.e., in near real-time) is available: Based on a model for stable historical behaviour
abnormal changes within newly acquired data can be detected. Different models are available
for modeling the stable historical behavior. A season-trend model (with harmonic seasonal pattern) is
  used as a default in the regresssion modelling.
}

\details{The package contains:
\itemize{
\item \code{\link[bfast]{bfast}}: Main function for iterative decomposition
  and break detection as described in Verbesselt et al (2010ab).
\item \code{\link[bfast]{bfastmonitor}}: Monitoring approach for detecting disturbances in
  near real-time (see Verbesselt et al. 2011, submitted to Remote Sensing and Environment).
\item \code{\link[bfast]{bfastpp}}: Data pre-processing for BFAST-type modeling.
\item Functions for plotting and printing, see  \code{\link[bfast]{bfast}}.
\item \code{\link[bfast]{simts}}: Artificial example data set.
\item \code{\link[bfast]{harvest}}: NDVI time series of a P. radiata plantation that is harvested.
\item \code{\link[bfast]{som}}: NDVI time series of locations in the south of Somalia to
  illustrate the near real-time disturbance approach
}
}

\author{
Jan Verbesselt [aut, cre], Achim Zeileis [aut], Rob Hyndman [ctb], Rogier De Jong [ctb]
}

\references{
Verbesselt J, Zeileis A, Herold M (2012).
  Near real-time disturbance detection using satellite image time series.
  \emph{Remote Sensing of Environment}, \bold{123}, 98--108.
  \url{http://dx.doi.org/10.1016/j.rse.2012.02.022}

Verbesselt J, Hyndman R, Newnham G, Culvenor D (2010).
  Detecting Trend and Seasonal Changes in Satellite Image Time Series.
  \emph{Remote Sensing of Environment}, \bold{114}(1), 106--115.
  \url{http://dx.doi.org/10.1016/j.rse.2009.08.014}

Verbesselt J, Hyndman R, Zeileis A, Culvenor D (2010).
  Phenological Change Detection while Accounting for Abrupt and Gradual Trends in Satellite Image Time Series.
  \emph{Remote Sensing of Environment}, \bold{114}(12), 2970--2980.
  \url{http://dx.doi.org/10.1016/j.rse.2010.08.003}

}

\keyword{ts}


text:
\name{print.comb_mat}
\alias{print.comb_mat}
\title{
Print the comb_mat Object
}
\description{
Print the comb_mat Object
}
\usage{
\method{print}{comb_mat}(x, ...)
}
\arguments{

  \item{x}{A combination matrix returned by \code{\link{make_comb_mat}}.}
  \item{...}{Other arguments}

}
\examples{
# There is no example
NULL

}

max_stars_repo_name: zhongmicai/complexHeatmap


text:
#Core level aggregation.
#Getting core, plot, site and global means and sd's for all observations made at the core level.
#This is currently soil pH, soil %C and soil C:N.
#clear environment, source paths.
rm(list=ls())
library(dplyr)
source('paths.r')
source('NEFI_functions/pC_uncertainty_neon.r')
source('NEFI_functions/cn_uncertainty_neon.r')
source('NEFI_functions/hierarch_core.means_JAGS.r')

#load data, subset and order.----
#samples with ITS data, along with DNA identifiers. Only take samples measured during peak greeness.
dp1.10801 <- readRDS(dp1.10108.00_output.path)
dp1.10801$geneticSampleID <- as.character(dp1.10801$geneticSampleID)
dp1.10801 <- dp1.10801[!(dp1.10801$geneticSampleID == ''),] #many samples missing geneticSampleID. throw them out.
dp1.10801 <- dp1.10801[order(dp1.10801$geneticSampleID),]

#moisture and pH data.
dp1.10086 <- readRDS(dp1.10086.00_output.path)
dp1.10086$site_date <- paste0(dp1.10086$siteID,'-',dp1.10086$dateID)
dp1.10086$geneticSampleID <- as.character(dp1.10086$geneticSampleID)

#soil C and N data.
dp1.10078 <- readRDS(dp1.10078.00_output.path)
dp1.10078$site_date_plot <- paste0(dp1.10078$siteID,'-',dp1.10078$dateID,'-',dp1.10078$plotID)
#remove analytical replicates
dp1.10078 <- dp1.10078[-which(dp1.10078$remarks %in% c("Duplicate samples", "Replicate samples", "Replicate samples.")),]
dp1.10078 <- dp1.10078[-which(dp1.10078$analyticalRepNumber == 2),]

#merge observations together.----
to_merge <- dp1.10086[,!c(colnames(dp1.10086) %in% colnames(dp1.10801))]
to_merge$geneticSampleID <- dp1.10086$geneticSampleID
merged <- merge(dp1.10801,to_merge)
to_merge <- dp1.10078[,!c(colnames(dp1.10078) %in% colnames(merged))]
to_merge$sampleID <- dp1.10078$sampleID

# fix chem data with N and C split up
to_merge$CNvals_merged <- FALSE
dup <- to_merge[duplicated(to_merge$sampleID),]
dup <- to_merge[to_merge$sampleID %in% dup$sampleID,]
dup$CNvals_merged <- TRUE
dup$testMethod <- "008 EA-IRMS organic d13C & C%, 009 EA-IRMS organic d15N & N%"
dup[,colnames(dup) %in% c("acidTreatment", "percentAccuracyQF", "instrument","analysisDate")]  <- NA # these rows lose their meaning once CN values are merged
shared_cols <- colnames(dup)[!colnames(dup) %in% c("organicCPercent", "nitrogenPercent", "CNratio")]
coalesce_by_column <- function(df) { # merge solution from: https://stackoverflow.com/questions/45515218/
  return(dplyr::coalesce(!!! as.list(df)))
}
reduced <- dup %>%
  group_by_at(vars(one_of(shared_cols))) %>%
  summarise_all(coalesce_by_column)
to_merge <- to_merge[!(to_merge$sampleID %in% reduced$sampleID),]
reduced <- reduced[,colnames(to_merge)] # reoder column names
reduced$CNratio <- round(reduced$organicCPercent/reduced$nitrogenPercent,1)
to_merge <- rbind(to_merge, as.data.frame(reduced)) # merge back together

#merge chem data with phys/DNA data
merged <- merge(merged,to_merge, by = 'sampleID', all.x=T)
merged$year <- substring(merged$dateID,1,4)
merged <- merged[!is.na(merged$siteID),]

#Subset to Peak Greenness and first observation. 581 observations. 13 sites.----
#FALSE - We subset to 2014 now, as Zoey realize this made things better.
#### THIS PEAK GREENNESS VS. 2014 THING MUST BE RESOLVED SOON.####
#But, that's not the first observation of everything- We need to choose molecular data separate?
#Determine which site year to grab. We want the first peak greenness observation for each site.
#merged <- merged[merged$sampleTiming == 'peakGreenness',]
merged <- merged[merged$year == '2014',]
merged$siteID <- as.character(merged$siteID)
sites <- unique(merged$siteID)
year_grab <- list()
for(i in 1:length(sites)){
  year_grab[[i]] <- min(merged[merged$siteID == sites[i], ]$year)
}
site_years <- data.frame(sites, unlist(year_grab))
colnames(site_years) <- c('siteID','year')
to_keep <- list()
for(i in 1:nrow(site_years)){
  to_keep[[i]] <- merged[merged$siteID == site_years$siteID[i] & merged$year == site_years$year[i],]
}
merged <- do.call(rbind, to_keep)
#get rid of duped geneticSampleIDs.
merged <- merged[!(duplicated(merged$geneticSampleID)),]

#finalize columns for core.level.----
#ATTN: WHY ISN'T ELEVATION IN THE ABOVE DATA PRODUCTS ANY LONGER?
#core.level <- merged[,c('sampleID','geneticSampleID','dnaSampleID','siteID','plotID','dateID','collectDate','horizon','elevation','soilMoisture','soilInWaterpH','organicCPercent','CNratio')]
core.level <- merged[,c('sampleID','geneticSampleID','dnaSampleID','siteID','plotID','dateID','collectDate','horizon',            'soilMoisture','soilInWaterpH','organicCPercent','CNratio')]
colnames(core.level)[(ncol(core.level) - 2) : ncol(core.level)] <- c('pH','pC','cn')
core.level$pC_sd <- pC_uncertainty_neon(core.level$pC)
core.level$cn_sd <- cn_uncertainty_neon(core.level$cn)

#get higher level observations and uncertainties (plot, site and global level.)----
#pC
pC.ag <- hierarch_core.means_JAGS(core.level$pC,core_plot = core.level$plotID)
pC.plot <- pC.ag$plot.table[,c('plotID','Mean','SD')]
pC.site <- pC.ag$site.table[,c('siteID','Mean','SD')]
pC.glob <- pC.ag$glob.table[,c('Mean','SD')]
colnames(pC.plot)[2:3] <- c('pC','pC_sd')
colnames(pC.site)[2:3] <- c('pC','pC_sd')
#C:N
cn.ag <- hierarch_core.means_JAGS(core.level$cn,core_plot = core.level$plotID)
cn.plot <- cn.ag$plot.table[,c('plotID','Mean','SD')]
cn.site <- cn.ag$site.table[,c('siteID','Mean','SD')]
cn.glob <- cn.ag$glob.table[,c('Mean','SD')]
colnames(cn.plot)[2:3] <- c('cn','cn_sd')
colnames(cn.site)[2:3] <- c('cn','cn_sd')
#pH
pH.ag <- hierarch_core.means_JAGS(core.level$pH,core_plot = core.level$plotID)
pH.plot <- pH.ag$plot.table[,c('plotID','Mean','SD')]
pH.site <- pH.ag$site.table[,c('siteID','Mean','SD')]
pH.glob <- pH.ag$glob.table[,c('Mean','SD')]
colnames(pH.plot)[2:3] <- c('pH','pH_sd')
colnames(pH.site)[2:3] <- c('pH','pH_sd')

#final output.
#core level.
core.out <- core.level

#y obs need to be in core.out. Essentially dropping any observation with no associated core-level covariate.
core.obs <- dp1.10801[dp1.10801$geneticSampleID %in% core.out$geneticSampleID,]
core.obs <- core.obs[!(duplicated(core.obs$geneticSampleID)),] #drop duplicated, as before.

#plot level.
plot.out <- merge(pC.plot , cn.plot, all = T)
plot.out <- merge(plot.out, pH.plot, all = T)
#Some plots in core-level not in plot-level.
to_add <- as.character(unique(core.out[!(core.out$plotID %in% plot.out$plotID),]$plotID))
to_add <- data.frame(to_add)
colnames(to_add) <- 'plotID'
plot.out <- plyr::rbind.fill(plot.out,to_add)
plot.out$siteID <- substring(plot.out$plotID, 1, 4)
plot.out <- plot.out[order(plot.out$plotID),]

#site level.
site.out <- merge(pC.site, cn.site, all = T)
site.out <- merge(site.out, pH.site, all = T)
#Some sites in plot-level not in site-level.
to_add <- as.character(unique(plot.out[!(plot.out$siteID %in% site.out$siteID),]$siteID))
to_add <- data.frame(to_add)
colnames(to_add) <- 'siteID'
site.out <- plyr::rbind.fill(site.out,to_add)

#global level
glob.out <- rbind(pC.glob, cn.glob, pH.glob)
pred <- c('pC','cn','pH')
glob.out <- cbind(pred,glob.out)

#save output.----
saveRDS(core.obs,  core_obs.path)
saveRDS(core.out, core_core.path)
saveRDS(plot.out, core_plot.path)
saveRDS(site.out, core_site.path)
saveRDS(glob.out, core_glob.path)


text:
#+---------------------------------------------------------------------------------------------------------------------+
#| 3. MAKE.LIB: (v. 2)                                                                                |
#|                                                                                |
#| This function makes the retention time library. It takes the following inputs in it's command:                      |
#| 1. The path to the edited evidence made by fix.evidence function.                                                   |
#| 2. The full name of the reference experiment as reported by find.ref.exp.                                           |
#| 3. The path for saving the library & Retention.time shifts info for experiments. We will need both of these files   |
#| for calculating the retention times of the peptide for each experiment.                                             |
#| The output is the Retention time library and the shifting coefficients for the experiments in the evidence file.    |
#| In the library, for each peptide mean, median, standard deviation and count of RTs is calculated.                   |
#|                                                                                |
#| This function needs 'MASS' package for the robust regression and it is assumed that it is installed and required.   |
#|                                                                                |
#| Example: make.lib(path.in = 'C:/../ev.ed.txt',                                                                      |
#|              path.out.lib = 'C:/../RT.lib.txt',                                                                     |
#|           path.out.coeffs = 'C:/../shift.coeffs.txt',                                                               |
#|                   ref.exp = '160623AEL160624PMRSAM00092_NC126')                                                     |
#| v. 2: Changes in how the function removes outliers.                                                                 |
#+---------------------------------------------------------------------------------------------------------------------+

make.lib <- function(path.in, path.out.lib, path.out.coeffs, ref.exp){
  #install.packages('MASS')
  #library('MASS')

  # load inputs
  ev.ed <- read.delim(path.in, header = TRUE, stringsAsFactors = FALSE)
  ref.pick <- ref.exp

  # load list of experiments
  exps <- unique(ev.ed$Raw.file)

  # load reference experiment data, and sort by peptide sequence
  reference.exp.l <- subset(ev.ed, ev.ed$PEP<5e-4 & ev.ed$Raw.file==ref.pick)
  reference.exp.l <- reference.exp.l[with(reference.exp.l, order(Sequence)), ]
  row.names(reference.exp.l) <- NULL

  # remove reference experiment from master list
  exps <- exps[!exps %in% ref.pick]

  # BUILDING THE LIBRARY:
  # initialize RT library data frame
  rt.lib <- data.frame(peptides=as.character(sort(unique(ev.ed$Sequence))),
                       matrix(NA, nrow=length(unique(ev.ed$Sequence)), ncol=length(unique(ev.ed$Raw.file))),
                       stringsAsFactors=FALSE)

  # load reference experiment but with less stringent threshold for PEP (0.05)
  reference.exp.h <- subset(ev.ed, ev.ed$PEP<0.05 & ev.ed$Raw.file==ref.pick)

  # load retention times in from reference experiment
  rt.lib$X1 <- reference.exp.h$Retention.time[match(rt.lib$peptides, reference.exp.h$Sequence)]

  # initialize shift data frame
  shift.coeffs <- data.frame(experiment=unique(ev.ed$Raw.file),
                             intercept=rep(NA, length(unique(ev.ed$Raw.file))),
                             coeff=rep(NA, length(unique(ev.ed$Raw.file))))
  # normalize to the reference experiment
  shift.coeffs[shift.coeffs$experiment==ref.pick, 2:3] <- 1

  # counters
  C10 <- 0
  C.res <- 0
  counter <- 2

  # for every experiment that is not the reference:
  for (i in exps) {
    counter <- counter+1

    cat('\r', 'Processing ', counter, '/', length(exps), ' ', i,  '                                          ')
    flush.console()

    # grab the experiment
    new.exp <- subset(ev.ed, ev.ed$Raw.file==i)
    row.names(new.exp) <- NULL

    # only take confidently identified peptides, sharing sequences w/ reference experiment, then sort
    new.exp.f <- subset(ev.ed, ev.ed$Raw.file==i & ev.ed$PEP<5e-4)
    new.exp.f <- new.exp.f[new.exp.f$Sequence %in% reference.exp.l$Sequence, ]
    new.exp.f <- new.exp.f[with(new.exp.f, order(Sequence)), ]
    row.names(new.exp.f) <- NULL
    new.exp.f <- droplevels(new.exp.f)

    # get data for these shared & confidently identified peptides
    ref.shared <- reference.exp.l[reference.exp.l$Sequence %in% new.exp.f$Sequence, ]
    row.names(ref.shared) <- NULL
    ref.shared <- droplevels(ref.shared)

    # iterate counter if less than 10 shared peptides
    # then continue loop to next experiment
    if (nrow(ref.shared)<10){
      C10 <- C10+1
      rt.lib[, counter] <- NA
      next
    }

    # build linear model and store coefficients
    fit <- rlm(ref.shared$Retention.time ~ new.exp.f$Retention.time, psi=psi.bisquare, maxit=100)
    shift.coeffs[shift.coeffs$experiment==i, 2:3] <- fit$coefficients

    fit.residuals <- as.matrix(data.frame(rep(1, nrow(new.exp.f)), new.exp.f$Retention.time)) %*% t(as.matrix(shift.coeffs[shift.coeffs$experiment==i, 2:3]))
    fit.residuals <- data.frame(peptides=new.exp.f$Sequence, shifted.RTs=as.numeric(fit.residuals))
    fit.residuals$ref.RTs <- ref.shared$Retention.time[match(fit.residuals$peptides, ref.shared$Sequence)]
    fit.residuals$residuals <- abs(fit.residuals$shifted.RTs-fit.residuals$ref.RTs)
    shift.coeffs[shift.coeffs$experiment==i, "mean.rt.res"] <- as.numeric(mean(fit.residuals$residuals))
    shift.coeffs[shift.coeffs$experiment==i, "median.rt.res"] <- as.numeric(median(fit.residuals$residuals))
    shift.coeffs[shift.coeffs$experiment==ref.pick, "mean.rt.res"] <- 0
    shift.coeffs[shift.coeffs$experiment==ref.pick, "median.rt.res"] <- 0

    edited.rts <- as.matrix(data.frame(rep(1, nrow(new.exp)), new.exp$Retention.time)) %*% t(as.matrix(shift.coeffs[shift.coeffs$experiment==i, 2:3]))
    edited.rts <- data.frame(peptides=new.exp$Sequence, shifted.RTs=as.numeric(edited.rts))
    #edited.rts$initial.RTs <- new.exp$Retention.time
    #edited.rts$residuals <- abs(edited.rts$shifted.RTs-edited.rts$initial.RTs)

    if (shift.coeffs[shift.coeffs$experiment==i, "median.rt.res"]>2){
      C.res <- C.res+1
      rt.lib[, counter] <- NA
      next
    }

    rt.lib[, counter] <- edited.rts$shifted.RTs[match(rt.lib$peptides, edited.rts$peptides)]
  }

  # LIBRARY POST PROCESSING:

  rt.lib[rt.lib<0] <- NA

  rt.lib$RT.count <- apply(rt.lib[,2:(ncol(rt.lib))], 1, function(x) sum(!is.na(x)))
  rt.lib <- rt.lib[!rt.lib$RT.count==0,]
  row.names(rt.lib) <- NULL

  rt.lib$rt.median <- apply(rt.lib[,2:(length(unique(ev.ed$Raw.file))+1)], 1, median, na.rm=T)

  rt.lib$rt.mean <- apply(rt.lib[,2:(length(unique(ev.ed$Raw.file))+1)], 1, mean, na.rm=T)

  rt.lib$sd.pep <- apply(rt.lib[,2:(length(unique(ev.ed$Raw.file))+1)], 1, sd, na.rm=T)
  rt.lib[is.na(rt.lib$sd.p), 'sd.pep'] <- 0

  write.table(rt.lib, path.out.lib, sep = '\t', row.names = FALSE, quote = FALSE)
  write.table(shift.coeffs, path.out.coeffs, sep = '\t', row.names = FALSE, quote = FALSE)
}


text:
# RICHIESTE?
prima_1 <- 3.178
prima_2 <- 6.39
seconda_1 <- 1.58
seconda_2 <- 3.178
terza <- 0.1011

regions <- c(-20, -1.615, -0.017, 1.58, 3.178, 4.775, 6.373, 7.97, 20)
f_values <- c(0, 0.1011, 0.3832, 0.7112, 0.9148, 0.9858, 0.999, 1)

# PROGRAMMA PRINCIPALE
#
#
plot.new()
EPSILON <- 0.001
x <- seq(-3, 9, by = EPSILON)

groups <- cut(x, regions)
fx <- f_values[groups]

par(new=TRUE)
plot(fx, x = x, type="l", xlab = "x", ylab = "F(x)", cex = EPSILON, col = "Blue4", bg = "transparent")

for(i in 2:length(regions)) {
  abline(v = regions[i], lty="dotted", col = "Gray")
  # abline(h = f_values[i], lty="dashed", col = "azure2")
  points(regions[i], f_values[i], col = "Blue4", bg = "Blue4", pch = 21, cex = 0.8)
  points(regions[i], f_values[i-1], col = "Blue4", pch = 21, cex = 0.8)
}

f <- function(x) {
  ifelse(x < regions[2], f_values[1],
    ifelse(x >= regions[2] && x < regions[3], f_values[2],
      ifelse(x >= regions[3] && x < regions[4], f_values[3],
        ifelse(x >= regions[4] && x < regions[5], f_values[4],
          ifelse(x >= regions[5] && x < regions[6], f_values[5],
            ifelse(x >= regions[6] && x < regions[7], f_values[6],
              ifelse(x >= regions[7] && x < regions[8], f_values[7],
                ifelse(x >= regions[8], f_values[8], NaN))))))))
}

# f_minus <- function(x) {
#   ifelse(x <= regions[2], f_values[1],
#     ifelse(x > regions[2] && x <= regions[3], f_values[2],
#       ifelse(x > regions[3] && x <= regions[4], f_values[3],
#         ifelse(x > regions[4] && x <= regions[5], f_values[4],
#           ifelse(x > regions[5] && x <= regions[6], f_values[5],
#             ifelse(x > regions[6] && x <= regions[7], f_values[6],
#               ifelse(x > regions[7] && x <=  regions[8], f_values[7],
#                  ifelse(x > regions[8], f_values[8], NaN))))))))
# }

f_minus <- function(x) {
  f(x - EPSILON)
}

print(f(prima_2)-f(prima_1))
print(f_minus(seconda_2)-f_minus(seconda_1))

for(i in 2:length(f_values)) {
  if(abs((f_values[i]-f_values[i-1]) - terza) < EPSILON)
    print(regions[i])
}

max_stars_repo_path: code/Esercizio 18.r

max_stars_repo_name: mfranzil/PSUniTN

max_issues_repo_path: code/Esercizio 18.r

max_issues_repo_name: mfranzil/PSUniTN

max_forks_repo_path: code/Esercizio 18.r

max_forks_repo_name: mfranzil/PSUniTN


text:
#' Clustring MCs based on the flow graph
#'
#' @param mct_id  mct network object
#'
#' @return A matrix on metacells and time with mean umi fraction over the gene module
#' @export
#'
mctnetwork_clust_flows = function(mct_id, K)
{
	mct = scdb_mctnetwork(mct_id)
	if(is.null(mct)) {
		stop("cannot find mctnet object ", mct_id, " when trying to plot net flows anchors")
	}
	mc = scdb_mc(mct@mc_id)

	fls = mctnetwork_get_flow_mat(mct,-1)
	fls3 = fls %*% fls %*% fls
	cr3 = tgs_cor(rbind(fls3, t(fls3)))
	diag(cr3) = 0
	cr3_2 = tgs_cor(cr3)
	diag(cr3_2) = 0

	hc = hclust(tgs_dist(t(cr3_2)), "ward.D")

	clf = cutree(hc, K)
	return(list(cmat = cr3_2, clust = clf, hc=hc))
}

mctnetwork_plot_cormat = function(mct_id, mat_id, gset_id, K, fn, w=1600, h=1600, text_cex=1, cls_pre_order=NULL, marks=NULL)
{
	mct = scdb_mctnetwork(mct_id)
	if(is.null(mct)) {
		stop("cannot find mctnet object ", mct_id, " when trying to plot net flows anchors")
	}
	clst_flows = mctnetwork_clust_flows(mct_id, K)

	cmat = clst_flows$cmat
	fclst = clst_flows$clust

	mc_ord = clst_flows$hc$order
	if(!is.null(cls_pre_order)) {
		cls_rank = 1:K
		names(cls_rank) = as.character(1:K)
		cls_rank[as.character(cls_pre_order)] = 1:K
		mc_mean_age = apply(mct@mc_t, 1, function(x) {return(mean(x*(1:length(x))/sum(x))) })
		mc_ord = order(-cls_rank[fclst]*1000+mc_mean_age)
	}

	mc = scdb_mc(mct@mc_id)

	shades = colorRampPalette(c("darkblue", "blue","white", "red", "yellow"))(1000)
	png(fn, w=w, h=h)
	n_mc = nrow(cmat)

	layout(matrix(c(1,2),nrow=2),heights=c(nrow(cmat)*9+50, 300))
	par(mar=c(0,5,4,5))
	image(cmat[mc_ord, mc_ord], zlim=c(-1,1), col=shades, yaxt='n', xaxt='n')
	N = length(mc_ord)

	mc_x = 1:length(mc_ord)
	names(mc_x) = 1:N
	mc_x[mc_ord] = 1:N
	for(i in 1:K) {
		abline(h=max(-0.5+mc_x[fclst == i])/(N-1))
		abline(v=max(-0.5+mc_x[fclst == i])/(N-1))
	}
	cl_x = tapply(mc_x, fclst, mean)
	cl_max = tapply(mc_x, fclst, max)

	mtext(1:K, side = 3, at=cl_x/N, las=1, cex=1.5)

	mtext((1:length(mc_ord))[mc_ord], side = 2, at=seq(0,1,l=n_mc), las=2, line = 2, cex=text_cex)
	mtext((1:length(mc_ord))[mc_ord], side = 4, at=seq(0,1,l=n_mc), las=2, line = 2, cex=text_cex)
	par(mar=c(3,5,0,5))
	image(as.matrix(mc_ord,nrow=1), col=mc@colors, yaxt='n', xaxt='n')
	mtext((1:length(mc_ord))[mc_ord], side = 1, at=seq(0,1,l=n_mc), las=2, line = 2, cex=text_cex)
	dev.off()

	if(!is.null(marks)) {
		egc = log2(mc@e_gc+1e-5)
		dir = sub(".png","_marks",fn)
		if(!dir.exists(dir)) {
			dir.create(dir)
		}
		for(g in marks) {
			png(sprintf("%s/%s.png", dir, g),w=1000,h=300)
			plot(1:length(mc_ord), egc[g, mc_ord], pch=19, col=mc@colors[mc_ord], ylab=g,xaxt='n')
			mtext(1:K, at=cl_x,side=1, las=2)

			abline(v=cl_max+0.5)
			grid()
			dev.off()
		}
	}
	mark_fn = sub("png", "mat.png", fn)
	mcell_mc_plot_marks(mct@mc_id, gset_id=gset_id, mat_id = mat_id, fig_fn=mark_fn, mc_ord = mc_ord, plot_cells=T)
	return(mc_ord)
}

max_stars_repo_name: ofirr/metacell

max_issues_repo_name: lbgbox/metacell

max_forks_repo_name: lbgbox/metacell


text:
library(data.table)
source("../summarise_phenotypes.r")

# Biomarkers for first 100K run
# n_chunks <- 1
# filename_root <- "biomarkers/pharma_exomes_biomarkers_parsed_output_100k_chunk."
# variable_info <- "variable-info/outcome_info_final_round3.tsv"
# coding_info_file <- "variable-info/data-coding-ordinal-info.txt"

# Note: coding info file for the first 100K run was the following file, rather than the recently updated file referred to below.
# coding_info_file <- "variable-info/data-coding-ordinal-info.txt"

n_exomes <- '200k'
n_chunks <- 10

variable_info <- "../variable-info/outcome_info_final_pharma_nov2019.tsv"
filename_root <- paste0("../../pharma_exomes_parsed_output_", n_exomes, "_chunk.")
coding_info_file <- "../variable-info/data-coding-ordinal-info-nov2019-update.txt"

# Read in all the data table codings
to_read <- paste("../WAS/codings/", dir("../WAS/codings"), sep="")
codings_tables <- list()

for(i in to_read) {
	name <- gsub("../WAS/codings/coding(.*).tsv", "\\1", i)
	codings_tables[[name]] <- read.table(i, header=TRUE, sep='\t', quote="", stringsAsFactors=FALSE)
}

for (i in 1:n_chunks)
{
	# Both sexes
	hist_filename <- paste0(filename_root, i, "_hist")
	pheno_summary <- paste0(filename_root, i, "_phenosummary.tsv")

	filename <- paste0(filename_root, i)
	tsv_filename <- paste0(filename, ".tsv")
	log_filename <- paste0(filename_root, i)
	log_file <- paste0(log_filename, ".log")
	tsv_data <- read.table(tsv_filename, header=TRUE, sep='\t')
	names(tsv_data)[1] <- "userId"

	outcome_info <- read.table(variable_info, sep='\t', quote="", comment.char="", header=TRUE)
	print(paste0("both sexes ", i))
	summary_file <-  get_hists_and_notes(hist_filename, tsv_data, log_file, outcome_info, codings_tables, samples_for_inclusion=TRUE, check=FALSE, start_column=4)
	write.table(summary_file, file=pheno_summary, col.names=TRUE, row.names=TRUE, sep='\t', quote=FALSE)
}

max_stars_repo_name: lganel/PHESANT

max_issues_repo_name: lganel/PHESANT

max_forks_repo_name: lganel/PHESANT


text:
#' Calculates the quantitative connectance, and the effective number of flows using the interaction matrix and a vector of species abundances
#'
#' The Quantitative connectance (Cq) takes into account both the distribution of per capita interaction strengths
#' among species in the web and the distribution of species’ abundances and quantifies the diversity of network fluxes,
#' if all the species have the same flux is equal to the directed connectance.
#' The mean or effective number of flows impinging upon or emanating from a tipical node (LDq) is based the average flow diversity, and when all flows are equal is similar to linkage density. Both measures are based in Shannon information theory.
#' The total interaction flux is measured
#' as ```T[i,j] <- d[i] * d[j] * interM[i,j]```. The effective Cq is calculated following the formulas in appendix 2 of [1], LDq follows [2]
#'
#' @references
#' 1. Fahimipour, A.K. & Hein, A.M. (2014). The dynamics of assembling food webs. Ecol. Lett., 17, 606–613
#'
#' 1. Ulanowicz, R.E. & Wolff, W.F. (1991). Ecosystem flow networks: Loaded dice? Math. Biosci., 103, 45–68
#'
#' @param interM per capita interaction strength matrix
#' @param d      species' abundances vector
#'
#' @return A list with Cq,the quantitative connectance index and LDq, =
#'
#' @aliases calcQuantitativeConnectance
#'
#' @export
#'
#' @examples
#'
#' # 3 predators 2 preys unequal fluxes
#' #
#' m <- matrix()
#' matrix(0,nrow=4,ncol=4)
#' m[1,2] <- m[1,3] <- m[3,4]<- .2
#' m[2,1] <- m[3,1] <- m[4,3] <- -2
#'
#' calc_quantitative_connectance(m, c(1,1,1,1))
#'
#'# Equal input and output fluxes
#'
#'m <- matrix(0,nrow=4,ncol=4)
#'m[1,2] <- m[1,3] <- m[3,4]<- 2
#'m[2,1] <- m[3,1] <- m[4,3] <- -2
#'calc_quantitative_connectance(m, c(1,1,1,1))

calc_quantitative_connectance <- function(interM,d){

  nsp <- length(d) # number of species
  if(ncol(interM)!=nrow(interM)) stop("interM has to be square, number of rows has to be equal to number of columns")
  if(nsp!=nrow(interM)) stop("length of species vector has to be equal to interM dimensions")

  TT <- matrix(0,nsp,nsp)
  SS <- matrix(0,nsp,nsp)

  Hout <- numeric(nsp)
  Hin  <- numeric(nsp)

  for(i in seq_along(d)) {
    for(j in seq_along(d)) {
      TT[i,j] <- abs(interM[i,j])*d[i]*d[j]
    }
  }
  TTcol <- colSums(TT)    # Input
  TTrow <- rowSums(TT)    # Output
  TTT <- sum(TTrow)

  for(k in seq_along(d)) {
    for(i in seq_along(d)) {
      SS[i,k] <- TT[i,k]/TTT * log2(TT[i,k]^2/(TTrow[i]*TTcol[k]))
      if(TT[i,k]!=0){
        Hin[k] <- Hin[k] - TT[i,k]/TTcol[k]*log(TT[i,k]/TTcol[k])
      }
      if(TT[k,i]!=0){
        Hout[k] <- Hout[k] - TT[k,i]/TTrow[k]*log(TT[k,i]/TTrow[k])
      }
    }
  }
  # SS1 is the \Sigma separated in input and output components see Ulanowicz eq 3
  # (SS1 <- (sum(TTcol*Hin)+sum(TTrow*Hout))/TTT) == sum(SS)
  #
  nout <- exp(Hout)
  nout[TTrow==0] <- 0
  nin <- exp(Hout)
  nin[TTcol==0] <- 0
  LDq <- 2^(-sum(SS,na.rm = T)/2)   # C= linkage Density / S
  Cq <- 1/(2*nsp^2)*(sum(nout+nin))  #
  return(list(Cq=Cq,LDq=LDq))
}

calcQuantitativeConnectance <- function(interM,d){
  calc_quantitative_connectance(interM,d)}

#' Calc the Quasi Sign Stability measure for antagonistic (predator-prey) or mgraph networks with multiple interactions.
#'
#' The proportion of matrices that are locally stable, these matrices are created by sampling the values of the community matrix
#' (the Jacobian) from a uniform distribution, preserving the sign structure [1]. If the 'ig' parameter is
#' an `mgraph` network it needs to have been built with the order `c("Competitive", "Mutualistic", "Trophic")`
#' It also calculates the mean of the real part of the maximum eingenvalue, which is also a measure of stability [2].
#' It uses a uniform distribution between 0 and maximum values given by the parameters `negative`, `positive` and `selfDamping`,
#' corresponding to the sign of interactions and self-limitation effect [3,4].
#' If the edges of the networks have a weigth attribute and `istrength` parameter is true, weigth will be used as interaction strength, then the limits of the uniform distribution
#' will be `negative*-x`, `positive*x` where x is the value of the weigth for the edge.
#' If the values of these parameters are 0 then there is no interaction of that kind.
#'
#' @references
#'
#' 1. Allesina, S. & Pascual, M. (2008). Network structure, predator - Prey modules, and stability in large food webs.
#' Theor. Ecol., 1, 55–64.
#' 2. Grilli, J., Rogers, T. & Allesina, S. (2016). Modularity and stability in ecological communities. Nat. Commun., 7, 12031
#' 3. Monteiro, A.B. & Del Bianco Faria, L. (2017). Causal relationships between population stability and food-web topology. Functional Ecology, 31, 1294–1300.
#' 4. Borrelli, J. J. 2015. Selection against instability: stable subgraphs are most frequent in empirical food webs. - Oikos 124: 1583–1588.

#'
#' @param ig  igraph or a list of igraph networks or mgraph network
#' @param nsim number of simulations to calculate QSS
#' @param ncores number of cores to use in parallel comutation if 0 it uses sequential processing
#' @param negative the maximum magnitude of the negative interaction (the effect of the predator on the prey) must be <= 0
#' @param positive the maximum magnitude of the positive interaction (the effect of the prey on the predator) must be >= 0
#' @param selfDamping the maximum magnitude of the self-limitation (the effect of the species on itself) must be <= 0
#' @param istrength If TRUE takes the weigth attribute of the network as interaction strength.
#'
#' @return a data.frame with the QSS, and MEing, the mean of the real part of the maximum eingenvalue
#'
#' @export
#'
#' @import igraph
#' @importFrom future sequential multiprocess
#' @importFrom future.apply future_lapply
#' @importFrom dplyr bind_rows
#'
#' @examples
#' \dontrun{
#'
#' g <- netData[[2]]
#'
#' tp <- calc_QSS(g)
#'
#' # Read Multiplex network and calculate QSS
#'
#' fileName <- c(system.file("extdata",  package = "multiweb"))
#' dn <- list.files(fileName,pattern = "^Kefi2015.*\\.txt$")
#' gt <- readMultiplex(dn,types,"inst/extdata", skipColumn = 2)
#' calc_QSS(gt)
#'
#' }
calc_QSS <- function(ig,nsim=1000,ncores=0,negative=-10, positive=0.1, selfDamping=-1,istrength=FALSE) {
  if(inherits(ig,"igraph")) {
    ig <- list(ig)
  } else if(class(ig[[1]])!="igraph") {
    stop("parameter ig must be an igraph object")
  }
  if(negative>0 || positive<0 || selfDamping>0)
    stop("Parameters should be negative<=0, positive>=0, selfDamping<=0")

  if(ncores) {
    cn <- future::availableCores()
    if(ncores>cn)
      ncores <- cn
    future::plan(multiprocess, workers=ncores)
    on.exit(future::plan(sequential))
  } else {
    future::plan(sequential)
  }

  if(class(ig)!='mgraph') {
    df <-  lapply(ig, function(red)
    {
      lred <- fromIgraphToMgraph(list(make_empty_graph(n=vcount(red)),make_empty_graph(n=vcount(red)),red),
                                   c("empty","empty","Antagonistic"))
      mat <- toGLVadjMat(lred,c("empty","empty","Antagonistic"),istrength = istrength)   #

      df <- future_lapply(seq_len(nsim), function(i){
        ranmat <- ranUnif(mat,negative,positive,selfDamping)
        eigs <- maxRE(ranmat)
      },future.seed = TRUE)
      df <-   do.call(rbind,df)
      data.frame(QSS=sum(df<0)/nsim,MEing=mean(df))
    })
  } else {
    mat <- toGLVadjMat(ig,c("Competitive", "Mutualistic", "Trophic"),istrength = istrength)   #
    df <- future_lapply(seq_len(nsim), function(i){
      ranmat <- ranUnif(mat,negative,positive,selfDamping)
      eigs <- maxRE(ranmat)
    },future.seed = TRUE)
    df <-   do.call(rbind,df)
    df <- data.frame(QSS=sum(df<0)/nsim,MEing=mean(df))
  }

  bind_rows(df)
}

# Auxiliar function of calc_QSS, Calculates a random community mattrix with a fixed signed structure
#
ranUnif <- function(motmat, negative=-10,positive=0.1,selfDamping=-1){
  newmat <- apply(motmat, c(1,2), function(x){
    #if(x==1){runif(1, 0, positive)}else if(x==-1){runif(1, negative, 0)} else{0}
    if(x>0){runif(1, 0, x*positive)}else if(x<0){runif(1,-x*negative, 0)} else{0}
  })
  diag(newmat) <- runif(nrow(motmat), selfDamping, 0)
  return(newmat)
}

# Auxiliar function of calc_QSS, Compute the eigenvalues and return the largest real part.
#
maxRE <- function(rmat){
  lam.max <- max(Re(eigen(rmat)$values))
  return(lam.max)
}

#' Function to calculate weighted network indices
#'
#' The function calculates: weighted linkage density, connectance, generality and vulnerability and the SDs of the last two based on [1]
#' and the level of omnivory, mean and maximum trophic level based on [2] using [NetIndices::TrophInd()] function.
#' The igraph networks must have the weight attribute.
#'
#' @param ig An igraph object with the weight attribute representing fluxes
#' @param ncores number of cores used to compute in parallel, if 0 sequential processing is used.
#'
#' @references
#'
#' 1. Bersier, LF. et al. (2002). Quantitative Descriptors Of Food-web Matrices. Ecology, 83(9), 2394–2407.
#' 2. Kones, JK. et al. (2009). Are network indices robust indicators of food web functioning? A Monte Carlo approach.
#'    Ecological Modelling, 220, 370–382.
#'
#'
#' @return a data.frame with the following fields:
#'
#'  \item{LD:}{ linkage density}
#'  \item{Connectance:}{ directed Connectance}
#'  \item{TLmean:}{ mean trophic level}
#'  \item{TLmax:}{ maximum trophic level}
#'  \item{LOmnivory:} { Level of omnivory, quantiﬁes mean of the variety in trophic levels of the preys of a consumer}
#'  \item{Vulnerability:}{ mean of number of consumers per prey}
#'  \item{VulSD:}{ the standard deviation of normalized Vulnerability}
#'  \item{Generality:}{ mean number of prey per consumer}
#'  \item{GenSD:}{ the standard deviation of normalized Generality}
#'
#'
#' @export
#' @importFrom NetIndices TrophInd
#' @import igraph
#' @importFrom foreach foreach %dopar%
#' @importFrom doFuture registerDoFuture
#' @importFrom future sequential multiprocess
#'
#' @examples
#'
#' # Generate a test network
#'
#' g <-   graph_from_literal( 1 -+ 4 -+ 7,2 -+ 5 -+7, 3-+6-+7, 7-+7, 4+-3, simplify = FALSE)
#'
#' # Add weight
#'
#' E(g)$weight <- sample(c(.1,.2,.8,.9),gsize(g),replace=TRUE)
#'
#' calc_weighted_topological_indices(g)

calc_weighted_topological_indices<- function(ig,ncores=0){

  if(inherits(ig,"igraph")) {
    ig <- list(ig)
  } else if(class(ig[[1]])!="igraph") {
    stop("parameter ig must be an igraph object")
  }

  registerDoFuture()
  if(ncores) {
    cn <- future::availableCores()
    if(ncores>cn)
      ncores <- cn
    future::plan(multiprocess, workers=ncores)
    on.exit(future::plan(sequential))
  } else {
    future::plan(sequential)
  }

  df <-  foreach(g=ig,.combine='rbind',.inorder=FALSE,.packages=c('igraph','NetIndices')) %dopar% {

    W.net <- as_adjacency_matrix(g,sparse=FALSE,attr="weight")
    res<-c()
    # The flux matrix

    ### Taxon-specific Shannon indices of inflows
    # sum of k species inflows --> colsums
    sum.in<-apply(W.net, 2, sum)

    # Diversity of k species inflows
    H.in.mat<- t(t(W.net)/sum.in)*t(log(t(W.net)/sum.in)) #columns divided by the total col sum
    H.in.mat[!is.finite(H.in.mat)] <- 0 #converts NaN to 0's
    H.in<- apply(H.in.mat, 2, sum)*-1

    ### Taxon-specific Shannon indices of outflows
    # sum of k speies outflows --> rowsums
    sum.out<-apply(W.net, 1, sum)

    # Diversity of k species outflows
    H.out.mat<- (W.net/sum.out)*log(W.net/sum.out) #rows divided by the total row sum
    H.out.mat[!is.finite(H.out.mat)] <- 0 #converts NaN to 0's
    H.out<- apply(H.out.mat, 1, sum)*-1

    # Effective number of prey or resources = N(R,k)
    # The reciprocal of H(R,k) --> N (R,k) is the equivalent number of prey for species k
    # N.res<-exp(H.in)
    N.res<-ifelse(sum.in==0, H.in, exp(H.in))

    # Effective number of predators or consumers = N(C,k)
    # The reciprocal of H(C,k) --> N (C,k) is the equivalent number of predators for species k
    # N.con<-exp(H.out)
    N.con<-ifelse(sum.out==0, H.out, exp(H.out))

    ### Quantitative Weighted and unweighted link density and weighted connectance
    no.species<-ncol(W.net)

    # The weighted link density (LDw) is:
    # In the weighted version the effective number of predators for species i is weighted by i's
    # contribution to the total outflow
    # the same is the case for the inflows

    tot.mat<- sum(W.net)
    # LD.w <- (sum((sum.in/tot.mat)*N.res) + sum((sum.out/tot.mat)*N.con))/2
    # equivalent to next formula, but next one is closer to manuscript
    res$qLD.w <- 1/(2*tot.mat)*(sum(sum.in*N.res) + sum(sum.out*N.con))

    # Weighted connectance
    res$qC.w<- res$qLD.w/no.species

    # defintion according to Bersier et al. 2002 top species = [0.99, 1]
    # int.sp<- pos.ind[pos.ind>0&pos.ind<0.99]# intermediate are all species that are not basal nor top

    # weighted quantitative Generality
    res$qG.w<-sum(sum.in*N.res/sum(W.net))

    # weighted quantitative Vulnerability
    res$qV.w<-sum(sum.out*N.con/sum(W.net))

    # Standard deviation of stanardized unweighted quantitative Generality
    s<-no.species # length(pos.ind)

    # Standard deviation of stanardized weighted quantitative Generality
    res$qGsd.w<-sd(s*sum.in*N.res/sum(N.res*sum.in))

    # Standard deviation of stanardized weighted quantitative Vulnerability
    res$qVsd.w<-sd(s*sum.out*N.con/sum(N.con*sum.out))

    # Weigthed trophic levels
    #
    wtl <- TrophInd(W.net)

    wOmn <- mean(wtl$OI)        # Omnivory Level

    data.frame( wLD=res$qLD.w, wConnectance=res$qC.w, wTLmean=mean(wtl$TL),wTLmax=max(wtl$TL),wLOmnivory=wOmn,
                wVulnerability=res$qV.w,wVulSD=res$qVsd.w,wGenerality=res$qG.w,wGenSD=res$qGsd.w)
  }
  return(df)

}

max_stars_repo_name: lsaravia/EcoNetwork

max_issues_repo_name: lsaravia/EcoNetwork

max_forks_repo_name: lsaravia/EcoNetwork


text:
# WORKFLOW FOR RMD
#
# For polygons, the only thing to change is the polygon dataset (poly_sf)
# and attribute to plot ("NONE" if no attribute)
# The rest of the lines stay the same.

poly_sf <-  BlueWhale_ImpHab_sf #EBSA_sf #NBNW_ImpHab_sf  #leatherback_sf #ClippedCritHab_sf #rockweed_sf
attribute <- "Activity"           #"NONE"         #"NONE"            #"RWP"

# THIS LINE WOULD BE THE SAME FOR ALL POLYGONS
# intersect polyData with 3 polygons (region, map area, and study box)
# master_intersect function is stored in fn_intersect_operations.R
clipped <- master_intersect(poly_sf, mapDataList, getRegion=TRUE)
# outputs: clipped$regionData, clipped$mapData, and clipped$studyData

# PLOT AREA MAP (using clipped$mapData)
if (!is.null(clipped$mapData)){
  plot_rr_sf(areaMap, clipped$mapData, attribute)
}

# PLOT REGION MAP (if needed, using clipped$regionData)
if (!is.null(clipped$regionData)){
  plot_rr_sf(regionMap, clipped$regionData, attribute)
}

max_stars_repo_name: AtlanticR/shinySpatialApp-SandBox

max_issues_repo_name: AtlanticR/shinySpatialApp-SandBox

max_forks_repo_name: AtlanticR/shinySpatialApp-SandBox


text:
library("tidyverse")
library("edgeR")
library("limma")
library("psych")
library("gridExtra")
library("scales")

# read in the prepped data
paw_spc <- read_tsv("QC_check_669x19.txt")

# save accessions in vector and remove from data table
accession <- paw_spc$Accession
paw_spc <- select(paw_spc, -Accession)

# replace zeros with 0.25
paw_spc[paw_spc == 0] <- 0.25
head(paw_spc)
nrow(paw_spc)

# function for simple normalization
SL_Norm <- function(df, color = NULL, plot = TRUE) {
    # This makes each channel sum to the average grand total
        # df - data frame of TMT intensities
        # returns a new data frame with normalized values

    # compute scaling factors to make colsums match the average sum
    norm_facs <- mean(c(colSums(df))) / colSums(df)
    cat("SL Factors:\n", sprintf("%-5s -> %f\n", colnames(df), norm_facs))
    df_sl  <- sweep(df, 2, norm_facs, FUN = "*")

    # visualize results and return data frame
    if(plot == TRUE) {
        boxplot(log10(df_sl), col = color, notch = TRUE, main = "SL Normalized data")
    }
    df_sl
}

# normalize the data before plotting
color <- c(rep('blue', 8), rep('red', 11))
paw_sl <- SL_Norm(paw_spc, color)

# load the data into edgeR data structures
# group labels need to be factors
group <- c(rep("Control", 8), rep("lepto", 11))
y <- DGEList(counts = paw_spc, group = group, genes = accession)

# run the TMM normalization (and library size corrections)
y <- calcNormFactors(y)

# check normalizations
y$samples

apply_tmm_factors <- function(y, color = NULL, plot = TRUE) {
    # computes the tmm normalized data from the DGEList object
        # y - DGEList object
        # returns a dataframe with normalized intensities

    # compute grand total (library size) scalings
    lib_facs <- mean(y$samples$lib.size) / y$samples$lib.size

    # the TMM factors are library adjustment factors (so divide by them)
    norm_facs <- lib_facs / y$samples$norm.factors
    cat("Overall Factors (lib.size+TMM):\n", sprintf("%-5s -> %f\n",
                                                     colnames(y$counts), norm_facs))

    # compute the normalized data as a new data frame
    df_tmm <- as.data.frame(sweep(y$counts, 2, norm_facs, FUN = "*"))
    colnames(df_tmm) <- str_c(colnames(y$counts), "_tmm")

    # visualize results and return data frame
    if(plot == TRUE) {
        boxplot(log10(df_tmm), col = color, notch = TRUE, main = "TMM Normalized data")
    }
    df_tmm
}

# get normalized data and check boxplots
paw_spc_tmm <- apply_tmm_factors(y, color)

# check clustering
plotMDS(y, col = color, main = "Control vs Lepto")

CV <- function(df) {
    # Computes CVs of data frame rows
        # df - data frame,
        # returns vector of CVs (%)
    ave <- rowMeans(df)    # compute averages
    sd <- apply(df, 1, sd) # compute standard deviations
    cv <- 100 * sd / ave   # compute CVs in percent (last thing gets returned)
}

labeled_boxplot <- function(df, ylim, title) {
    # Makes a box plot with the median value labeled
        # df - data frame with data to compute CVs of
        # ylim - upper limit for y-axis
        # title - plot title
    cv = CV(df)
    boxplot(cv, ylim = c(0, ylim), notch = TRUE, main = title)
    text(x = 0.65, y = boxplot.stats(cv)$stats[3],
         labels = round(boxplot.stats(cv)$stats[3], 1))
}

# define the columns for each condition
C <- 1:8
L <- 9:19
CD <- c(1, 2, 5, 6, 7, 8)
LD <- c(9, 11, 12, 13, 15, 16, 17, 19)

# see what effect TMM had on CV distributions
par(mfrow = c(2, 2))
labeled_boxplot(paw_spc[C], 200, "All Control before")
labeled_boxplot(paw_spc[L], 200, "All Lepto before")
labeled_boxplot(paw_spc_tmm[C], 200, "All Control after")
labeled_boxplot(paw_spc_tmm[L], 200, "All Lepto after")

# dropping 2 controls and 3 lepto samples
par(mfrow = c(2, 2))
labeled_boxplot(paw_spc[CD], 200, "Control before")
labeled_boxplot(paw_spc[LD], 200, "Lepto before")
labeled_boxplot(paw_spc_tmm[CD], 200, "Control after")
labeled_boxplot(paw_spc_tmm[LD], 200, "Lepto after")
par(mfrow = c(1, 1))

# load the data into edgeR data structures
# group labels need to be factors
countsd <- paw_spc[c(C, LD)]
colord <- c(rep('blue', 8), rep('red', 8))
groupd <- c(rep("Control", 8), rep("lepto", 8))
yd <- DGEList(counts = countsd, group = groupd, genes = accession)

# run the TMM normalization (and library size corrections)
yd <- calcNormFactors(yd)

# check normalizations
yd$samples

# get normalized data and make boxplots
paw_spc_tmm <- apply_tmm_factors(yd, colord)

# check clustering
plotMDS(yd, col = colord, main = "ControlD vs LeptoD")

# log the session
sessionInfo()