使用作者程式碼重複結果

2020 年 3 月 30 日
筆記

序言

第三單元第十二+十三講：使用作者程式碼重複結果課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53

這一篇會是程式碼密集型，因為原文作者的程式碼真的寫的很長！

下載作者的Github

https://github.com/KPLab/SCS_CAF

文件布局如下：

下載好以後，需要將那兩個tar.gz文件解壓縮

看第二個R腳本 Processing.R

讀入表達量數據

# 首先指定操作路徑  Path_Main<-"~/scrna/SCS_CAF-master"  # 然後讀入原始的第一個細胞板數據  plate1_raw<-read.delim(paste(Path_Main,"/SS2_15_0048/counts.tab",sep=""),header=TRUE,check.names=FALSE,sep="t")  > plate1_raw[1:3,1:3]      gene A3 A6  1 Adora1  0  0  2  Sntg1  0  0  3  Prim2  0  0

作者這裡考慮到重複基因名的問題

# 的確存在重複基因名  > length(as.character(plate1_raw$gene))  [1] 24490  > length(unique(as.character(plate1_raw$gene)))  > sum(duplicated(as.character(plate1_raw$gene)))  [1] 492  [1] 23998  # 用這個查看  as.character(plate1_raw$gene)[duplicated(as.character(plate1_raw$gene))]    # 看一下make.unique的用法  > make.unique(c("a", "a"))  [1] "a"   "a.1"  # 將重複基因名變為唯一的名字  plate1_raw$gene<-make.unique(as.character(plate1_raw$gene))  > sum(duplicated(as.character(plate1_raw$gene)))  [1] 0  # 對樣本重新命名  colnames(plate1_raw)[2:length(colnames(plate1_raw))]<-paste("SS2_15_0048_",colnames(plate1_raw)[2:length(colnames(plate1_raw))],sep="")

同樣的，對於0049板，也是上述操作，最後將它們按照gene這一列進行合併，並把gene轉為行名

expr_raw<-merge(plate1_raw,plate2_raw,by="gene",all=TRUE)  rownames(expr_raw)<-as.character(expr_raw$gene)  expr_raw$gene<-NULL

最後計算一下dropout的比例（結果有點高）：

# 計算dropout的比例  sum(expr_raw==0)/(dim(expr_raw)[1]*dim(expr_raw)[2])  # 0.8305757

讀入ERCC數據

# 也是類似上面?的操作  # plate1  plate1_raw_ercc<-read.delim(paste(Path_Main,"/SS2_15_0048/counts-ercc.tab",sep=""),header=TRUE,check.names=FALSE,sep="t")  plate1_raw_ercc$gene<-make.unique(as.character(plate1_raw_ercc$gene))  colnames(plate1_raw_ercc)[2:length(colnames(plate1_raw_ercc))]<-paste("SS2_15_0048_",colnames(plate1_raw_ercc)[2:length(colnames(plate1_raw_ercc))],sep="")    # plate2  plate2_raw_ercc<-read.delim(paste(Path_Main,"/SS2_15_0049/counts-ercc.tab",sep=""),header=TRUE,check.names=FALSE,sep="t")  plate2_raw_ercc$gene<-make.unique(as.character(plate2_raw_ercc$gene))  colnames(plate2_raw_ercc)[2:length(colnames(plate2_raw_ercc))]<-paste("SS2_15_0049_",colnames(plate2_raw_ercc)[2:length(colnames(plate2_raw_ercc))],sep="")    # 最後合併、計算ERCC dropout  expr_raw_ercc<-merge(plate1_raw_ercc,plate2_raw_ercc,by="gene",all=TRUE)  rownames(expr_raw_ercc)<-as.character(expr_raw_ercc$gene)  expr_raw_ercc$gene<-NULL    sum(expr_raw_ercc==0)/(dim(expr_raw_ercc)[1]*dim(expr_raw_ercc)[2])  # 0.6267691

看一下ERCC在各個細胞的表達量分布：

barplot(sort(as.numeric(colSums(expr_raw_ercc)),decreasing=TRUE),ylab="SPIKE LIBRARY SIZE",xlab="CELL INDEX")

然後做一個直方圖，把一定數量的樣本中ERCC表達量合併作一個bin：

hist(log2(as.numeric(colSums(expr_raw_ercc))+1),col="brown",       main="Distribution of Spike Library Sizes",       xlab="Spike Library Size",breaks=20)

將內源基因與ERCC spike-in合併

先看看分別有多少：

> print(paste0("There are ",nrow(expr_raw)," endogenous genes"))  [1] "There are 24490 endogenous genes"  > print(paste0("There are ",nrow(expr_raw_ercc)," spikes"))  [1] "There are 92 spikes"

合併起來：

all.counts.raw<-rbind(expr_raw,expr_raw_ercc)  > dim(all.counts.raw)  [1] 24582   768

然後重新計算dropout的比例：

sum(all.counts.raw==0)/(dim(all.counts.raw)[1]*dim(all.counts.raw)[2])  # 0.8298129

一共有7153個基因在所有細胞中表達量均為0：

dim(all.counts.raw[rowSums(all.counts.raw)==0,])  # 7153  768

關於原文去掉的52個細胞

根據一些指標去掉了52個細胞

作者也把這52個細胞的質控結果讀入了R：

cell_QC<-read.delim(paste(Path_Main,"/qc/qc_2plates.filtered_cells.txt",sep=""),row.names=1,header=TRUE,sep="t")  > dim(cell_QC)  [1] 52  6

在原始矩陣中也要去掉這些細胞：

rownames(cell_QC)<-gsub("__","_",rownames(cell_QC))    all.counts.raw<-subset(all.counts.raw,select=colnames(all.counts.raw)[!colnames(all.counts.raw)%in%rownames(cell_QC)])  > dim(all.counts.raw)  [1] 24582   716

過濾細胞後，重新拆分成count矩陣和ERCC矩陣：

# 得到原始count矩陣  expr_raw<-subset(expr_raw,select=colnames(expr_raw)[!colnames(expr_raw)%in%rownames(cell_QC)])  # 得到ERCC矩陣  expr_raw_ercc<-subset(expr_raw_ercc,select=colnames(expr_raw_ercc)[!colnames(expr_raw_ercc)%in%rownames(cell_QC)])

分別對count矩陣和ERCC矩陣過濾

all.counts.raw<-all.counts.raw[rowMeans(all.counts.raw)>0,]  expr_raw<-expr_raw[rowMeans(expr_raw)>=1,]  # count矩陣過濾後只剩下10835個基因  expr_raw_ercc<-expr_raw_ercc[rowMeans(expr_raw_ercc)>0,]  # ERCC也有原來的92個變成了89個

然後畫CV vs Mean圖

library("matrixStats")  # 首先還是計算CV值  mean_expr_raw<-as.numeric(rowMeans(expr_raw,na.rm=TRUE))  sd_expr_raw<-rowSds(as.matrix(expr_raw),na.rm=TRUE)  cv_squared_expr_raw<-(sd_expr_raw/mean_expr_raw)^2  # plot函數中使用(縱坐標~橫坐標)  plot(log10(cv_squared_expr_raw)~log10(mean_expr_raw),       pch=20,cex=0.5,xlab="log10 ( mean raw count )",       ylab="log10 ( CV^2)",main="RAW COUNTS")  # 接下來添加ERCC的資訊(畫上紅點)  mean_expr_raw_ercc<-as.numeric(rowMeans(expr_raw_ercc,na.rm=TRUE))  sd_expr_raw_ercc<-rowSds(as.matrix(expr_raw_ercc),na.rm=TRUE)  cv_squared_expr_raw_ercc<-(sd_expr_raw_ercc/mean_expr_raw_ercc)^2  points(log10(cv_squared_expr_raw_ercc)~log10(mean_expr_raw_ercc),col="red",pch=20,cex=1.5)  # 然後對ERCC添加loess擬合曲線  fit_expr_raw_ercc<-loess(log10(cv_squared_expr_raw_ercc)[is.finite(log10(mean_expr_raw_ercc))]~log10(mean_expr_raw_ercc)[is.finite(log10(mean_expr_raw_ercc))],span=1)  # 從小到大排個序  j<-order(log10(mean_expr_raw_ercc)[is.finite(log10(mean_expr_raw_ercc))])  #  lines(fit_expr_raw_ercc$fitted[j]~log10(mean_expr_raw_ercc)[is.finite(log10(mean_expr_raw_ercc))][j],col="red",lwd=3)

又根據擬合結果，進行了預測得到期望值，然後過濾得到符合期望CV值的基因，最後只留下5316個基因：

pred_expr_raw<-predict(fit_expr_raw_ercc,log10(mean_expr_raw))  filtered_expr_raw<-expr_raw[log10(cv_squared_expr_raw)>=pred_expr_raw,]  filtered_expr_raw<-filtered_expr_raw[grepl("NA",rownames(filtered_expr_raw))==FALSE,]  > dim(filtered_expr_raw)  [1] 5316  716

可以看到，它的過濾從原來24490的基因，然後過濾掉沒表達的基因剩10835個，然後又需要符合期望，剩5000多個。最後就是拿這5000多個基因做下游分析

看第三個R腳本 Dimensionality_reduction.R

這個腳本需要RPKM結果，因此需要先跑完上面第二個完整的腳本

降維主要使用tSNE，聚類使用dbscan（它的作用和hclust或者kmeans差不多）

上來先跑50次tSNE：

library(Rtsne)  N_tsne <- 50  tsne_out <- list(length = N_tsne)  KL <- vector(length = N_tsne)  set.seed(1234)  for(k in 1:N_tsne)  {    tsne_out[[k]]<-Rtsne(t(log10(RPKM+1)),initial_dims=30,verbose=FALSE,check_duplicates=FALSE,                         perplexity=27, dims=2,max_iter=5000)    KL[k]<-tail(tsne_out[[k]]$itercosts,1)    print(paste0("FINISHED ",k," TSNE ITERATION"))  }  names(KL) <- c(1:N_tsne)  # 可以看到這裡選擇最小的KL作為50次中效果最優的tSNE，然後主要關注tsne結果的itercosts  opt_tsne <- tsne_out[[as.numeric(names(KL)[KL==min(KL)])]]$Y  opt_tsne_full<-tsne_out[[as.numeric(names(KL)[KL==min(KL)])]]  save(tsne_out,opt_tsne,opt_tsne_full,file="step3-tsne-out.Rdata")

然後使用dbscan聚類：

library(dbscan)  plot(opt_tsne,  col=dbscan(opt_tsne,eps=3.1)$cluster, pch=19, xlab="tSNE dim 1", ylab="tSNE dim 2")

如果使用kmeans方法：

plot(opt_tsne,  col=kmeans(opt_tsne,centers = 4)$clust, pch=19, xlab="tSNE dim 1", ylab="tSNE dim 2")

看看這兩種聚類方法的相關性:

> table(kmeans(opt_tsne,centers = 4)$clust,dbscan(opt_tsne,eps=3.1)$cluster)          0   1   2   3   4    1   1 226   0   0   0    2   0   0 144   0   0    3   0  27   0   0  44    4   0 231   0  43   0  # 左側是kmeans，上方是dbscan。

發現有一個點是離群值，所以把它放到細胞數量最多的那個組：

library(dbscan)  CAFgroups<-dbscan(opt_tsne,eps=3.1)$cluster  CAFgroups_full<-dbscan(opt_tsne,eps=3.1)  CAFgroups[CAFgroups==0]<-1  CAFgroups_full$cluster[CAFgroups_full$cluster==0]<-1  plot(opt_tsne,  col=CAFgroups, pch=19, xlab="tSNE dim 1", ylab="tSNE dim 2")

(補充)了tSNE，還可以對PCA可視化：

CAFgroups<-dbscan(opt_tsne,eps=3.1)$cluster  CAFgroups_full<-dbscan(opt_tsne,eps=3.1)  CAFgroups[CAFgroups==0]<-1  CAFgroups_full$cluster[CAFgroups_full$cluster==0]<-1    RPKM.PCA<-prcomp(log2(t(RPKM)+1), center=TRUE)  plot(RPKM.PCA$x,main="first PCA", pch=19, col=CAFgroups)

其實有了上面這個tSNE聚類圖，我們就能把基因的表達量映射上去，很像Seurat的FeaturePlot()做的那樣。但是這裡作者自己創造函數（參考第五個腳本：Plotting.R）

需要用到基因名、表達量矩陣、tsne坐標

plot.feature2<-function(gene, tsne.output=tsne.out, DATAuse=DATA){    plot.frame<-data.frame(x=tsne.output$Y[,1], y=tsne.output$Y[,2], log2expr=as.numeric(log2(DATAuse[gene,]+1)))    p<-ggplot(plot.frame,aes(x=x, y=y, col=log2expr))+      geom_point(size=1) +      labs(title=paste(gene))+      theme_classic()+      scale_color_gradientn(colors = c("#FFFF00", "#FFD000","#FF0000","#360101"), limits=c(0,14))+      theme(axis.title = element_blank())+      theme(axis.text = element_blank())+      theme(axis.line = element_blank())+      theme(axis.ticks = element_blank())+      theme(plot.title = element_text(size=20,face="italic"))+      theme(legend.title  = element_blank())+      theme(legend.position = "none")    return(p)  }    library(ggplot2)  opt_tsne <- tsne_out[[as.numeric(names(KL)[KL==min(KL)])]]$Y  opt_tsne_full<-tsne_out[[as.numeric(names(KL)[KL==min(KL)])]]  load(file='RPKM.full.Rdata')  load(file='CAFgroups.Rdata')  plot.feature2("Pdgfra", opt_tsne_full, RPKM.full)

另外小提琴圖的程式碼更是長：它是用來繪製不同基因的表達量在不同聚類分組的差異

需要用到基因名、表達量矩陣、tsne坐標

plot.violin2 <- function(gene, DATAuse, tsne.popus, axis=FALSE, legend_position="none", gene_name=FALSE){    testframe<-data.frame(expression=as.numeric(DATAuse[paste(gene),]), Population=tsne.popus$cluster)    testframe$Population <- as.factor(testframe$Population)    colnames(testframe)<-c("expression", "Population")      col.mean<-vector()    for(i in levels(testframe$Population)){      col.mean<-c(col.mean,mean(testframe$expression[which(testframe$Population ==i)]))    }    col.mean<-log2(col.mean+1)    col.means<-vector()    for(i in testframe$Population){      col.means<-c(col.means,col.mean[as.numeric(i)])    }    testframe$Mean<-col.means    testframe$expression<-log2(testframe$expression+1)      p <- ggplot(testframe, aes(x=Population, y=expression, fill= Mean, color=Mean))+      geom_violin(scale="width") +      labs(title=paste(gene), y ="log2(expression)", x="Population")+      theme_classic() +        scale_color_gradientn(colors = c("#FFFF00", "#FFD000","#FF0000","#360101"), limits=c(0,14))+      scale_fill_gradientn(colors = c("#FFFF00", "#FFD000","#FF0000","#360101"), limits=c(0,14))+      theme(axis.title.y =  element_blank())+      theme(axis.ticks.y =  element_blank())+      theme(axis.line.y =   element_blank())+      theme(axis.text.y =   element_blank())+      theme(axis.title.x = element_blank())+      theme(legend.position=legend_position )      if(axis == FALSE){      p<-p+        theme( axis.line.x=element_blank(),               axis.text.x = element_blank(),               axis.ticks.x = element_blank())    }    if(gene_name == FALSE){      p<-p+  theme(plot.title = element_blank())    }else{ p<-p + theme(plot.title = element_text(size=10,face="bold"))}    p  }    # 例如  plot.violin2(gene = "Pdgfra", DATAuse = RPKM.full, tsne.popus = CAFgroups_full)

看第四個R腳本 Differential_gene_expression.R

主要利用了ROTS包（Reproducibility-optimized test statistic），對每個亞群和其他幾個亞群共同體進行比較

差異分析重點就在：表達矩陣和分組資訊

library(ROTS)  library(plyr)  # 首先針對第一亞群和其他亞群比較(把其他亞群定義為234)  groups<-CAFgroups  groups[groups!=1]<-234    ROTS_input<-RPKM.full[rowMeans(RPKM.full)>=1,]  ROTS_input<-as.matrix(log2(ROTS_input+1))  # 運行程式碼很簡單，重點就是data和group參數  results_pop1 = ROTS(data = ROTS_input, groups = groups , B = 1000 , K = 500 , seed = 1234)  # 最後根據FDR值得到第一組和其他組比較的差異基因  summary_pop1<-data.frame(summary(results_pop1, fdr=1))  head(summary_pop1)  ##         Row ROTS.statistic pvalue FDR  ## Rgs5   8345      -19.89479      0   0  ## Higd1b 4559      -17.49991      0   0  ## Abcc9   393      -16.44638      0   0  ## Pdpn   7193       16.02262      0   0  ## Fbln2  3635       15.80534      0   0  ## Rgs4   8344      -15.62123      0   0    # 同理，對第2組可以與1、3、4合併組比較；對第3組可以和第1、2、4組比較；對第4組可以和第1、2、3組比較  # 都得到以後，共同保存  save(summary_pop1,summary_pop2,summary_pop3,summary_pop4,       file = 'ROTS_summary_pop.Rdata')

每個亞群可以挑top18基因繪製熱圖

population_subset<-c(rownames(summary_pop1[summary_pop1$ROTS.statistic<0,])[1:18],rownames(summary_pop2[summary_pop2$ROTS.statistic<0,])[1:18],rownames(summary_pop3[summary_pop3$ROTS.statistic<0,])[1:18],rownames(summary_pop4[summary_pop4$ROTS.statistic<0,])[1:18])  RPKM_heatmap<-RPKM.full[population_subset,]    RPKM_heatmap<-RPKM_heatmap[,order(CAFgroups_full$cluster)]  RPKM_heatmap<-log2(RPKM_heatmap+1)    popul.col<-sort(CAFgroups_full$cluster)  popul.col<-replace(popul.col, popul.col==1,"#1C86EE" )  popul.col<-replace(popul.col, popul.col==2,"#00EE00" )  popul.col<-replace(popul.col, popul.col==3,"#FF9912" )  popul.col<-replace(popul.col, popul.col==4,"#FF3E96" )  library(gplots)    #pdf("heatmap_genes_population.pdf")  heatmap.2(as.matrix(RPKM_heatmap),ColSideColors = as.character(popul.col), tracecol = NA, dendrogram = "none",col=bluered, labCol = FALSE, scale="none", key = TRUE, symkey = F, symm=F,  key.xlab = "", key.ylab = "", density.info = "density", key.title = "log2(RPKM+1)", keysize = 1.2, denscol="black", Colv=FALSE)

當然，原文還使用了其他幾種差異分析方法，放在這裡，可以做日後參考

################################  ####### 第一種：DESeq2 ###########  ################################  library("scran")  library("limSolve")  library(scater)  library(DESeq2)  ann<-data.frame(Plate = factor(unlist(lapply(strsplit(colnames(RPKM.full),"_"),function(x) x[3]))), Population = factor(gsub("(3|4)","2",as.character(CAFgroups)),levels=c("1","2")))  ann<-data.frame(Population = factor(gsub("(3|4)","2",as.character(CAFgroups)),levels=c("1","2")))  rownames(ann)<-colnames(RPKM.full)    ddsFullCountTable <- DESeqDataSetFromMatrix(    countData = all.counts.raw[rownames(RPKM.full),],    colData = ann,    design = ~ Population)    ddsFullCountTable<-DESeq(ddsFullCountTable)  DESeq_result<-results(ddsFullCountTable)  DESeq_result<-DESeq_result[order(DESeq_result$padj, DESeq_result$pvalue),]  head(DESeq_result,30)  write.table(DESeq_result[grep("ERCC", rownames(DESeq_result), invert=TRUE),], "DESeq_result.txt", col.names = TRUE, row.names = TRUE, quote = FALSE, sep="t")    ################################  ####### 第二種：EdgeR ############  ################################  library("edgeR")  edgeR_Data<-DGEList(counts=all.counts.raw[rownames(RPKM.full),], group=ann$Population)  edgeR_Data<-estimateCommonDisp(edgeR_Data)  edgeR_Data<-estimateTagwiseDisp(edgeR_Data)  edgeR_result<-exactTest(edgeR_Data)  edgeR_result_table<-edgeR_result$table  edgeR_result_table<-edgeR_result_table[order(edgeR_result_table$PValue),]  edgeR_result_table<-edgeR_result_table[grep("ERCC", rownames(edgeR_result_table), invert=TRUE),]  write.table(edgeR_result_table, "EdgeR_result.txt", col.names = TRUE, row.names = TRUE, quote=FALSE, sep="t")    ################################  ###### 第三種：Wilcox ##########  ################################  NumPerm<-1000  POP1_expr<-subset(RPKM.full,select=rownames(ann)[ann==1])  POP2_expr<-subset(RPKM.full,select=rownames(ann)[ann==2])  p_wilcox<-vector()  p_t<-vector()  p_perm<-vector()  statistics<-vector()  median_POP1_expr<-vector()  median_POP2_expr<-vector()  a<-seq(from=0,to=length(rownames(RPKM.full)),by=1000)    print("START DIFFERENTIAL GENE EXPRESSION BETWEEN POP1 AND POP2")  for(i in 1:length(rownames(RPKM.full)))  {    p_wilcox<-append(p_wilcox,wilcox.test(as.numeric(POP1_expr[rownames(RPKM.full)[i],]),as.numeric(POP2_expr[rownames(RPKM.full)[i],]))$p.value)    statistics<-append(statistics,wilcox.test(as.numeric(POP1_expr[rownames(RPKM.full)[i],]),as.numeric(POP2_expr[rownames(RPKM.full)[i],]))$statistic)    p_t<-append(p_t,t.test(as.numeric(POP1_expr[rownames(RPKM.full)[i],]),as.numeric(POP2_expr[rownames(RPKM.full)[i],]))$p.value)    p_perm<-append(p_perm,PermTest_Median(as.numeric(POP1_expr[rownames(RPKM.full)[i],]),as.numeric(POP2_expr[rownames(RPKM.full)[i],]),NumPerm))    median_POP1_expr<-append(median_POP1_expr,median(as.numeric(POP1_expr[rownames(RPKM.full)[i],])))    median_POP2_expr<-append(median_POP2_expr,median(as.numeric(POP2_expr[rownames(RPKM.full)[i],])))    if(i%in%a){print(paste("FINISHED ",i," GENES",sep=""))}  }  fold_change<-median_POP1_expr/median_POP2_expr  log2_fold_change<-log2(fold_change)  p_adj<-p.adjust(p_wilcox,method="fdr")  output_wilcox<-data.frame(GENE=rownames(RPKM.full),POP1_EXPR=median_POP1_expr,POP234_EXPR=median_POP2_expr,FOLD_CHANGE=fold_change,LOG2FC=log2_fold_change,WILCOX_STAT=statistics,P_T_TEST=p_t,P_PERM=p_perm,P_WILCOX=p_wilcox,FDR=p_adj)  output_wilcox<-output_wilcox[order(output_wilcox$P_PERM,output_wilcox$FDR,output_wilcox$P_WILCOX,output_wilcox$P_T_TEST,-abs(output_wilcox$LOG2FC)),]  print(head(output_wilcox,20))  write.table(output_wilcox,file="Wilcox_Perm_de_results.txt",col.names=TRUE,row.names=FALSE,quote=FALSE,sep="t")    ################################  ####### 第四種：SCDE WORKFLOW  ################################  library("scde")  # factor determining cell types  sg<-factor(as.numeric(CAFgroups))  # the group factor should be named accordingly  names(sg)<-colnames(expr_raw)  table(sg)  # define two groups of cells  #groups<-sg  groups <- factor(gsub("(3|4)","2",as.character(sg)),levels=c("1","2"))  table(groups)  # calculate models    cd<-apply(expr_raw,2,function(x) {storage.mode(x) <- 'integer'; x})  colnames(cd)<-colnames(expr_raw)  o.ifm<-scde.error.models(counts=cd,groups=groups,n.cores=4,threshold.segmentation=TRUE,save.crossfit.plots=FALSE,save.model.plots=FALSE,verbose=1)  print(head(o.ifm))  # filter out cells that don't show positive correlation with  # the expected expression magnitudes (very poor fits)  valid.cells<-o.ifm$corr.a > 0  table(valid.cells)  o.ifm<-o.ifm[valid.cells, ]  # estimate gene expression prior  o.prior<-scde.expression.prior(models=o.ifm,counts=cd,length.out=400,show.plot=FALSE)  # run differential expression tests on all genes.  ediff<-scde.expression.difference(o.ifm,cd,o.prior,groups=groups,n.randomizations=100,n.cores=4,verbose=1) #batch=batch  # top upregulated genes  ediff_order<-ediff[order(abs(ediff$Z),decreasing=TRUE), ]  head(ediff_order,20)  write.table(ediff_order,file="scde_de_results_1_vs_234.txt",col.names=TRUE,row.names=TRUE,quote=FALSE,sep="t")

Tags: