Script-代码分析--FPKM

count2FPKM

• commandArgs，通过命令行传递参数。将以上你的脚本保存为test.R文件，然后再终端运行Rscript test.R args1 args2 args3, args1-3分别是要输入的参数。

• quit() 退出程序

1. 参数设置

# commandArgs用来接收命令行参数，并赋值给args
args<-commandArgs(T)
# 对传递进来的参数进行条件判断，满足条件则继续执行脚本，不满足则退出，并打印出正确调用的提示。
if(length(args)<2){
cat("Usages Rscript count2FPKM.r  gene_count.tsv gene_length.txt gene_FPKM.tsv","\n")
quit('no')
}

2. 脚本函数语句

# count --> FPKM的函数
# 输入：count matrix; gene长度文件
# 输出：FPKM matrix
calc_fpkm <- function(counts, gene2length) {
    ## counts:gene-sample matrix; gene2length:gene-genelength txt file
    ## match() 返回输入的counts矩阵，gene2length文件中，能匹配上的gene2length的index，并用order保存index 文件
    order <- match(rownames(counts),gene2length[,1])
    ## 筛选与gene2length对应的counts矩阵
    norm <- counts[!is.na(order),]
    order <- order[!is.na(order)]
    ## 计算FPKM
    ## 基因长度与测序深度的加权
    weighted_sum_lengths <- colSums(norm*as.numeric(gene2length[order,2]))
    fpkms <- t(t(norm)/weighted_sum_lengths)*10^9
    ## 保留4位小数
    fpkms <- round(fpkms,4)
    return(fpkms)
}

或者

calc_fpkm <- function(counts, gene2length){
    genes <- gene2length[,1]  %in% rownames(counts)
    sub_counts <- counts[genes,]
    sub_gene2length<- gene2length[genes,]
    weighted_sum_lengths <- colSums(sub_counts*as.numeric(sub_gene2length[genes,2]))
    fpkms <- t(t(sub_counts )/weighted_sum_lengths)*10^9
    fpkms <- round(fpkms,4)
    return(fpkms)
}

3. 脚本输出结果保存
   函数调用

counts<-read.table(args[1],head=T,sep="\t",row.names=1,check.names=F)
gene2length<-read.table(args[2],sep="\t",head=T,check.names=F)
res<-calc_fpkm(counts,gene2length)

结果保存

## 筛选在细胞中有表达的genes
res<-res[rowSums(res)>0,]
## 转化为dataframe格式
### 一般 data.frame（data, check.names=F）和 write.table(data,quote=F,sep="\t",row.names=F) 一起使用。
res <- data.frame(gene=rownames(res),res,check.names=F)
write.table(res,args[3],quote=F,sep="\t",row.names=F)

脚本中R语言的函数/方法

• match() 数据筛选
  返回x中元素在table中的位置, 返回index
  match(x, table, nomatch = NA_integer_, incomparables = NULL)
  返回x中每个元素在table中是否存在, 返回TRUR/FALSE
  x %in% table
• check.names = FALSE
  这个坑经常会遇到，一般对于sampleID barcode等，会在index中出现“ .”, 在write.table()/read.table()/fread() 等读取或者保存数据文件的时，如果不设置check.names=FALSE, 就会将“.”, 改成“-”， 从而影响后续的数据分析。
• rowSums
  这个函数在数据筛选中经常用到，count[rowSums(count)>0,] 用于筛选至少在一个样本/细胞中有表达的gene。