Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

change import function in R #39

Open
wants to merge 17 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
planemo/
15 changes: 7 additions & 8 deletions tools/regionalgam/ab-index.R
Original file line number Diff line number Diff line change
@@ -1,14 +1,13 @@
#!/usr/bin/env Rscript

args = commandArgs(trailingOnly=TRUE)
args <- commandArgs(trailingOnly = TRUE)
source(args[1])

input <- data.frame(data.table::fread(args[2]))
pheno <- read.table(args[3], header = TRUE, sep = " ")

tryCatch({input = read.table(args[2], header=TRUE,sep=" ")},finally={input = read.table(args[2], header=TRUE,sep=",")})
pheno = read.table(args[3], header=TRUE,sep=" ")

if("TREND" %in% colnames(input)){
input <- input[input$TREND==1,c("SPECIES","SITE","YEAR","MONTH","DAY","COUNT")]
if ("TREND" %in% colnames(input)) {
input <- input[input$TREND == 1, c("SPECIES", "SITE", "YEAR", "MONTH", "DAY", "COUNT")]
}
data.index <- abundance_index(input, pheno)
write.table(data.index, file="data.index", row.names=FALSE, sep=" ")
data_index <- abundance_index(input, pheno)
write.table(data_index, file = "data_index", row.names = FALSE, sep = " ")
12 changes: 7 additions & 5 deletions tools/regionalgam/ab-index.xml
Original file line number Diff line number Diff line change
Expand Up @@ -3,20 +3,22 @@
<macros>
<import>regionalgam_macros.xml</import>
</macros>
<expand macro="rg_r_requirements"/>
<requirements>
<requirement type="package" version="1.12.2">r-data.table</requirement>
</requirements>
<command detect_errors="exit_code"><![CDATA[
Rscript '$__tool_directory__/ab-index.R'
'$__tool_directory__/dennis-gam-initial-functions.R'
'$count_file' '$input2'
'$output'
'$__tool_directory__/dennis-gam-initial-functions.R'
'$count_file' '$input2'
'$output'
]]>
</command>
<inputs>
<expand macro="rg_count_file"/>
<param format="tabular" name="input2" type="data" label="Regional expected pattern of abundance" help="Flight curve output, tabular file."/>
</inputs>
<outputs>
<data format="tabular" name="output" from_work_dir="data.index" />
<data format="tabular" name="output" from_work_dir="data_index" />
</outputs>
<tests>
<test>
Expand Down
10 changes: 5 additions & 5 deletions tools/regionalgam/autocorr-res-acf.R
Original file line number Diff line number Diff line change
@@ -1,10 +1,10 @@
#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)
args <- commandArgs(trailingOnly = TRUE)
load(args[1])

png('output-acf.png',width = 480, height = 480, units = "px")
graph<-acf(residuals(mod,type="normalized"),plot=FALSE)
plot(graph, main="Acf residuals")
png("output_acf.png", width = 480, height = 480, units = "px")
graph <- acf(residuals(mod, type = "normalized"), plot = FALSE)
plot(graph, main = "Acf residuals")
invisible(dev.off())
sink("output-acf.txt")
sink("output_acf.txt")
print(graph)
12 changes: 6 additions & 6 deletions tools/regionalgam/autocorr-res-acf.xml
Original file line number Diff line number Diff line change
Expand Up @@ -6,17 +6,17 @@
<expand macro="rg_r_requirements"/>
<command detect_errors="exit_code"><![CDATA[
Rscript '$__tool_directory__/autocorr-res-acf.R'
'$gls_model'
'$output'
'$output_res_values'
'$gls_model'
'$output'
'$output_res_values'
]]>
</command>
<inputs>
<expand macro="rg_gls_model"/>
</inputs>
<outputs>
<data format="png" name="output" from_work_dir="output-acf.png" />
<data format="txt" name="output_res_values" from_work_dir="output-acf.txt" />
<data format="png" name="output" from_work_dir="output_acf.png" />
<data format="txt" name="output_res_values" from_work_dir="output_acf.txt" />
</outputs>
<tests>
<test>
Expand All @@ -31,7 +31,7 @@
</edam_topics>
<help><![CDATA[
======================================
Model residuals autocorrelation check
Model residuals autocorrelation check
======================================

This tool is an implementation of the autocorr-res_acf function from RegionalGAM package: https://github.com/RetoSchmucki/regionalGAM/
Expand Down
182 changes: 91 additions & 91 deletions tools/regionalgam/dennis-gam-initial-functions.R

Large diffs are not rendered by default.

8 changes: 4 additions & 4 deletions tools/regionalgam/flight-curve.R
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)
args <- commandArgs(trailingOnly = TRUE)
source(args[1]) #TODO replace by library(regionalGAM) if available as official package from bioconda
tryCatch({input = read.table(args[2], header=TRUE,sep=" ")},finally={input = read.table(args[2], header=TRUE,sep=",")})
dataset1 <- input[,c("SPECIES", "SITE", "YEAR", "MONTH", "DAY", "COUNT")]
input <- data.frame(data.table::fread(args[2]))
dataset1 <- input[, c("SPECIES", "SITE", "YEAR", "MONTH", "DAY", "COUNT")]
pheno <- flight_curve(dataset1, MinVisit = args[3], MinOccur = args[4])

write.table(pheno, file="pheno", row.names=FALSE, sep=" ")
write.table(pheno, file = "pheno", row.names = FALSE, sep = " ")
15 changes: 8 additions & 7 deletions tools/regionalgam/flight-curve.xml
Original file line number Diff line number Diff line change
Expand Up @@ -5,14 +5,15 @@
</macros>
<requirements>
<requirement type="package" version="1.8_28">r-mgcv</requirement>
<requirement type="package" version="1.12.2">r-data.table</requirement>
</requirements>
<command detect_errors="exit_code"><![CDATA[
Rscript '$__tool_directory__/flight-curve.R'
'$__tool_directory__/dennis-gam-initial-functions.R'
'$count_file'
$minvisit
$minoccur
'$output'
'$__tool_directory__/dennis-gam-initial-functions.R'
'$count_file'
$minvisit
$minoccur
'$output'
]]>
</command>
<inputs>
Expand All @@ -34,7 +35,7 @@
<has_n_columns n="6"/>
<has_text_matching expression="&quot;Pyronia&#032;tithonus&quot;&#009;2003&#009;[0-9]&#009;[0-9]&#009;[0-9]&#009;[0-9]"/>
</assert_contents>
</output>
</output>
</test>
</tests>
<edam_topics>
Expand All @@ -50,7 +51,7 @@ This tool is an implementation of the flight_curve function RegionalGAM package:

This function computes the annual phenology on a specific region.

The output tabular file can be used to impute expected count values with the abundance index computation tool.
The output tabular file can be used to impute expected count values with the abundance index computation tool.

|

Expand Down
14 changes: 7 additions & 7 deletions tools/regionalgam/glmmpql.R
Original file line number Diff line number Diff line change
Expand Up @@ -2,15 +2,15 @@
library(nlme)
library(MASS)

args = commandArgs(trailingOnly=TRUE)
input = read.table(args[1], header=TRUE,sep=" ")
glmm.mod_fullyear <- glmmPQL(regional_gam~ as.factor(YEAR)-1,data=input,family=quasipoisson,random=~1|SITE, correlation = corAR1(form = ~ YEAR | SITE),verbose = FALSE)
args <- commandArgs(trailingOnly = TRUE)
input <- read.table(args[1], header = TRUE, sep = " ")
glmm_mod_fullyear <- glmmPQL(regional_gam ~ as.factor(YEAR) - 1, data = input, family = quasipoisson, random = ~ 1 | SITE, correlation = corAR1(form = ~ YEAR | SITE), verbose = FALSE)

col.index <- as.numeric(glmm.mod_fullyear$coefficients$fixed)
col_index <- as.numeric(glmm_mod_fullyear$coefficients$fixed)
year <- unique(input$YEAR)

write.table(col.index, file="output-glmmpql", row.names=FALSE, sep=" ")
write.table(col_index, file = "output-glmmpql", row.names = FALSE, sep = " ")

png('output-plot.png')
plot(year,col.index,type='o', xlab="year",ylab="collated index")
png("output-plot.png")
plot(year, col_index, type = "o", xlab = "year", ylab = "collated index")
invisible(dev.off())
14 changes: 7 additions & 7 deletions tools/regionalgam/gls-adjusted.R
Original file line number Diff line number Diff line change
Expand Up @@ -2,15 +2,15 @@
library(nlme)
library(MASS)

args = commandArgs(trailingOnly=TRUE)
input1 = read.table(args[1], header=TRUE) #input1=col.index =glmmpql-output
input2 = read.table(args[2], header=TRUE,sep=" ") #input2=data.index =abundance_index-output
args <- commandArgs(trailingOnly = TRUE)
input1 <- read.table(args[1], header = TRUE)
input2 <- read.table(args[2], header = TRUE, sep = " ")

input1<-as.matrix(input1)
input1 <- as.matrix(input1)

year <- unique(input2$YEAR)
mod <- gls(input1~year, correlation = corARMA(p=2))
summary<-summary(mod)
mod <- gls(input1 ~ year, correlation = corARMA(p = 2))
summary <- summary(mod)

save(mod, file = "mod_adjusted.rda")
capture.output(summary, file="mod_adjusted-summary.txt")
capture.output(summary, file = "mod_adjusted-summary.txt")
14 changes: 7 additions & 7 deletions tools/regionalgam/gls.R
Original file line number Diff line number Diff line change
Expand Up @@ -2,15 +2,15 @@
library(nlme)
library(MASS)

args = commandArgs(trailingOnly=TRUE)
input1 = read.table(args[1], header=TRUE) #input1=col.index =glmmpql-output
input2 = read.table(args[2], header=TRUE,sep=" ") #input2=data.index =abundance_index-output
args <- commandArgs(trailingOnly = TRUE)
input1 <- read.table(args[1], header = TRUE)
input2 <- read.table(args[2], header = TRUE, sep = " ")

input1<-as.matrix(input1)
input1 <- as.matrix(input1)

year <- unique(input2$YEAR)
mod <- gls(input1~year)
summary<-summary(mod)
mod <- gls(input1 ~ year)
summary <- summary(mod)

save(mod, file = "mod.rda")
capture.output(summary, file="mod-summary.txt")
capture.output(summary, file = "mod-summary.txt")
20 changes: 10 additions & 10 deletions tools/regionalgam/plot-trend.R
Original file line number Diff line number Diff line change
Expand Up @@ -2,22 +2,22 @@
library(nlme)
library(MASS)

args = commandArgs(trailingOnly=TRUE)
input = read.table(args[1], header=TRUE,sep=" ") #input=data.index =ab_index-output
args <- commandArgs(trailingOnly = TRUE)
input <- read.table(args[1], header = TRUE, sep = " ")
load(args[2])

glmm.mod_fullyear <- glmmPQL(regional_gam~ as.factor(YEAR)-1,data=input,family=quasipoisson,random=~1|SITE, correlation = corAR1(form = ~ YEAR | SITE),verbose = FALSE)
glmm_mod_fullyear <- glmmPQL(regional_gam ~ as.factor(YEAR) - 1, data = input, family = quasipoisson, random = ~ 1 | SITE, correlation = corAR1(form = ~ YEAR | SITE), verbose = FALSE)

col.index <- as.numeric(glmm.mod_fullyear$coefficients$fixed)
col_index <- as.numeric(glmm_mod_fullyear$coefficients$fixed)
year <- unique(input$YEAR)

png('output-plot-trend.png')
plot(year,col.index, type='o',xlab="year",ylab="collated index")
abline(mod,lty=2,col="red")
png("output-plot-trend.png")
plot(year, col_index, type = "o", xlab = "year", ylab = "collated index")
abline(mod, lty = 2, col = "red")
invisible(dev.off())

d<-data.frame(year,col.index)
colnames(d)<-c("year","collated index")
write.table(d,"s_trend.tsv",sep=" ",row.names=FALSE)
d <- data.frame(year, col_index)
colnames(d) <- c("year", "collated index")
write.table(d, "s_trend.tsv", sep = " ", row.names = FALSE)

cat("Show trend line on abundance plot")