R/SVM.R

In s18692001/rSeqTU: rSeqTU package

#' Using Support Vector Machine to train and generate TU prediction results
#' @import Rsubread Rsamtools QuasR e1071 seqinr grid gridBase ggplot2 reshape2
#' @import plyr caret mlbench Gviz GenomicRanges
#' @param positive_training The file of positive training dataset.
#' @param negative_training The file of negative training dataset.
#' @param positive_strand_testing The file of positive strand testing dataset.
#' @param negative_strand_testing The file of negative strand testing dataset.
#' @param file_RNAseqSignals The .NA file generated by previous step
#' @param file_gff The .gff file of reference genome
#' @param output_prefix The prefix of output file name
#' @param genome_name The file of referecne genome
#' @return TU prediction results
#' @export


# Import training and testing dataset for SVM
TU_SVM <- function(positive_training, negative_training, positive_strand_testing, negative_strand_testing, file_RNAseqSignals,
                   file_gff, output_prefix, genome_name){

positive_training <- read.delim("SimulatedPositiveTUMatrix.txt")
negative_training <- read.delim("SimulatedNegativeTUMatrix.txt")
positive_strand_testing <- read.delim("TargetPositiveTUMatrix.txt")
negative_strand_testing <- read.delim("TargetNegativeTUMatrix.txt")

pos_train_df <- as.data.frame(positive_training)
str(pos_train_df)
neg_train_df <- as.data.frame(negative_training)
str(neg_train_df)
pos_str_test_df <- as.data.frame(positive_strand_testing)
neg_str_test_df <- as.data.frame(negative_strand_testing)


# add labels in the last column, 1 is presence and 0 is absence of TU
pos_train_df$target_variable <- rep(1,nrow(pos_train_df))
neg_train_df$target_variable <- rep(0,nrow(neg_train_df))
train_data <- rbind(pos_train_df, neg_train_df)

# Convert target variable to factors
train_data[["target_variable"]] <- factor(train_data[["target_variable"]])

# Remove rows with NA
train_data <- train_data[complete.cases(train_data),]

# Dataset summarized details
summary(train_data)

# Slice data into training and testing dataset
intrain <- createDataPartition(y = train_data$target_variable, p= 0.7, list = FALSE)
training <- train_data[intrain,]
testing <- train_data[-intrain,]


### Feature Selection

# prepare training scheme
control <- trainControl(method="repeatedcv", number=10, repeats=3)
# train the model
model <- train(target_variable~., data=train_data, method="lvq", preProcess="scale", trControl=control)
# estimate variable importance
importance <- varImp(model, scale=FALSE)
# summarize importance
print(importance)
# plot importance
plot(importance)
imp_df = varImp(model)$importance
selected_feature = rownames(imp_df[order(imp_df$X0, decreasing = TRUE),])[1:8]
factors = c(selected_feature)

# define the control using a random forest selection function
control <- rfeControl(functions=rfFuncs, method="cv", number=10)
# run the RFE algorithm
results <- rfe(train_data[,1:17], train_data[,18], sizes=c(1:17), rfeControl=control)
# summarize the results
print(results)
# list the chosen features
predictors(results)
# plot the results
plot(results, type=c("g", "o"))


featurePlot(x = train_data[, 1:17],
            y = train_data$target_variable ,
            plot = "box",
            strip=strip.custom(par.strip.text=list(cex=.7)),
            scales = list(x = list(relation="free"),
                          y = list(relation="free")))

featurePlot(x = train_data[, 1:17],
            y = train_data$target_variable,
            plot = "density",
            strip=strip.custom(par.strip.text=list(cex=.7)),
            scales = list(x = list(relation="free"),
                          y = list(relation="free")))

f1 <- "target_variable"

# Training the Linear SVM model
trctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 3)
svm_Linear <- train(as.formula(paste(f1, paste(factors, collapse="+"), sep = "~")), method = "svmLinear", data = training,
                    trControl=trctrl,
                    preProcess = c("center", "scale"),
                    tuneLength = 10)
svm_Linear


# Testing set prediction
test_pred <- predict(svm_Linear, newdata = testing)
test_pred

# Confusion matrix
confusionMatrix(test_pred, testing$target_variable)

# Use to model to predict on real dataset
#real_predict <- predict(svm_Linear, newdata = real_data)
#real_predict
pos_str_test_df <- replace(pos_str_test_df, is.na(pos_str_test_df), 0)
SVMForward.Target.scaled.Predict <- predict(svm_Linear, newdata = pos_str_test_df)
neg_str_test_df <- replace(neg_str_test_df, is.na(neg_str_test_df), 0)
SVMReverse.Target.scaled.Predict <- predict(svm_Linear, newdata = neg_str_test_df)

# Read .NA file
SelectedRNAseqData<-read.table(file_RNAseqSignals, head=F)

# Read gff file
x=read.delim(file_gff, skip=7, header = F)
x1<-x[x$V3=="gene",]
x2<-within(x1[1,], y<-data.frame(do.call('rbind',strsplit(as.character(x1[1,]$V9), ';', fixed=TRUE))))
x3<-within(x2, y<-data.frame(do.call('rbind',strsplit(as.character(x2$y$X2), '=', fixed=TRUE))))
Ngene = data.frame(x3$y$X2, x3$V4, x3$V5, x3$V7)
for (i in 2:nrow(x1)){
  x2<-within(x1[i,], y<-data.frame(do.call('rbind',strsplit(as.character(x1[i,]$V9), ';', fixed=TRUE))))
  x3<-within(x2, y<-data.frame(do.call('rbind',strsplit(as.character(x2$y$X2), '=', fixed=TRUE))))
  x4<-data.frame(x3$y$X2, x3$V4, x3$V5, x3$V7)
  Ngene = rbind(Ngene,x4)
}
colnames(Ngene)=c("ID", "Start", "End", "Strand") # Set Column Names
Ngene=Ngene[order(Ngene[,2]),] #Order the list by start positions of each gene

#Select Genes
SelectedForwardGenes=Ngene[Ngene[,4]=="+",] #Select all genes on forward strand
dim(SelectedForwardGenes)
SelectedReverseGenes=Ngene[Ngene[,4]=="-",] #Select all genes on reverse strand
dim(SelectedReverseGenes)
endIndexForwardTargetPredict=nrow(SelectedForwardGenes)-1

# Initialize the signal matrix for the forward Results
MatrixRow = length(SVMForward.Target.scaled.Predict)
MatrixColumn = 1
MatrixForward = mat.or.vec(MatrixRow, MatrixColumn+2)

# Assign values to the signal matrix
for (i in 1:MatrixColumn){
  MatrixForward[,i] = as.matrix(SVMForward.Target.scaled.Predict)
}
mode(MatrixForward) = 'numeric'
# parse the frequency
for (i in 1:MatrixRow){
  positiveNumber = length(MatrixForward[i,MatrixForward[i,]==1])
  negativeNumber = MatrixColumn - positiveNumber
  if (positiveNumber>=negativeNumber){
    MatrixForward[i,MatrixColumn+1] = 1
  }else{
    MatrixForward[i,MatrixColumn+1] = -1
  }
  MatrixForward[i,MatrixColumn+2] = max(positiveNumber,negativeNumber)/MatrixColumn
}

# initialize the signal matrix for the reverse Results
MatrixRow = length(SVMReverse.Target.scaled.Predict)
MatrixReverse = mat.or.vec(MatrixRow, MatrixColumn+2)

# assign values to the signal matrix
for (i in 1:MatrixColumn){
  MatrixReverse[,i] = as.matrix(SVMReverse.Target.scaled.Predict)
}
mode(MatrixReverse) = 'numeric'

# Parse the frequency
for (i in 1:MatrixRow){
  positiveNumber = length(MatrixReverse[i,MatrixReverse[i,]==1])
  negativeNumber = MatrixColumn - positiveNumber
  if (positiveNumber>=negativeNumber){
    MatrixReverse[i,MatrixColumn+1] = 1
  }else{
    MatrixReverse[i,MatrixColumn+1] = -1
  }
  MatrixReverse[i,MatrixColumn+2] = max(positiveNumber,negativeNumber)/MatrixColumn
}

ForwardTargetTUsPrediction <- as.matrix(MatrixForward[,MatrixColumn+1])
mode(ForwardTargetTUsPrediction) <- 'numeric'
ReverseTargetTUsPrediction <- as.matrix(MatrixReverse[,MatrixColumn+1])
mode(ReverseTargetTUsPrediction) <- 'numeric'
ForwardTargetTUsPredictionFreq <- as.matrix(MatrixForward[,MatrixColumn+2])
mode(ForwardTargetTUsPredictionFreq) <- 'numeric'
ReverseTargetTUsPredictionFreq <- as.matrix(MatrixReverse[,MatrixColumn+2])
mode(ReverseTargetTUsPredictionFreq) <- 'numeric'
ForwardTUtable<-cbind(GenerateFinalTUTable(ForwardTargetTUsPrediction,SelectedForwardGenes,ForwardTargetTUsPredictionFreq),Strand="+")
ReverseTUtable<-cbind(GenerateFinalTUTable(ReverseTargetTUsPrediction,SelectedReverseGenes,ReverseTargetTUsPredictionFreq),Strand="-")
head (ForwardTUtable)
head (ReverseTUtable)
AllTUtable<-merge(ForwardTUtable, ReverseTUtable, all=TRUE)
AllTUtable[,2] <- paste(AllTUtable[,1], AllTUtable[,2], sep="")
AllTUtable <- AllTUtable[,c(3,4,6,5,2)]
AllTUtable <- AllTUtable[order(AllTUtable$FinalTUOnePosStart),]
row.names(AllTUtable) <- paste("TU",formatC(c(1:nrow(AllTUtable)), width=4,flag="0"), sep="")
FinalTUTableFileName = paste(output_prefix, "_FinalTUTable_forPlot", sep="")
write.table(AllTUtable, file = FinalTUTableFileName, sep = "\t", col.names=FALSE)

system(paste("sed -e 's/\"//g' ",FinalTUTableFileName," > ",FinalTUTableFileName,".clean",sep=""))
system(paste("mv ",FinalTUTableFileName,".clean ",FinalTUTableFileName,sep=""))


## Conversion of result table to gff format
M9_result_table <- read.delim(FinalTUTableFileName, header = FALSE)
M9 <- as.data.frame(M9_result_table)
M9$val <- ifelse(M9$V4 == '+',yes = 1,no= -1)
M9$seqname <- rep(genome_name, nrow(M9))
head(M9)
M9 <- M9[,c(8,2,3,7)]
head(M9)
sink(paste(output_prefix, "_result.bedgraph", sep=""))
cat("track color=0,255,0 altColor=255,0,0\n")
sink()
write.table(M9, file=paste(output_prefix, "_result.bedgraph", sep=""), sep="\t", quote=FALSE, row.names = FALSE, col.names = FALSE, append=TRUE)

}