minor base提取代码

从全基因组minor alelle中提取基因区minor base

#此类的txt文件，若直接用read.csv导入，并不能进行行列操作(非[i,j])，因为此时的数据是由\t间隔的字符串；

#用strsplit("file","symble")命令将向量中的每个元素间隔开；

#用do.call(rbind/cbind,file)来将向量转换为行列的格式，便于之后的操作；

#用which(condition)能将文件中复合条件的行提取出来；

#在linux下对gff文件进行初步处理；(在操作数据文件前，先熟悉文件格式和各个变量的含义)

#提取minoralelle中的minorbase：

chr1=strsplit(chr,"\t")

chr2=do.call(rbind,chr1)

x=nrow(chr2)

for(i in 1:x){
        lin=chr2[i,]

        chrom=chr2[i,3]

        pos=chr2[i,4]

        C1=length(which((lin=="CC"))) #此四行能再进一步简化

        T1=length(which((lin=="TT")))

        A1=length(which((lin=="AA")))

        G1=length(which((lin=="GG")))

        a=c(C1,G1,T1,A1)

        b=which(a/500<=0.2&a/500!=0)

        if(length(b)!=0){

        cat(chrom,pos,b,"\n",append=TRUE,file="improvesnp.txt")

}     

}

(#之后用R处理时，改进用cat(content,"\n",append=TRUE,file="")命令来将每次循环的结果分别输出，"\n"换行符，append=TRUE将输出内容接在末尾)

(#通过cat的改进，使得代码运算时间相较之前减少了10X，且能即时查看当前循环输出文件的结果)

#将提取出来的minorbase与gff基因注释文件的位置进行比较：

chr1=do.call(rbind,strsplit(readLines("egchr1.txt")," "))

chr2=read.csv("pos1.txt",stringsAsFactors = F)

x=nrow(chr1)

y=nrow(chr2)

for(i in 1:x){

lin1=as.numeric(chr1[i,2])

for(j in 1:y){

if(lin1>=chr2[j,2] & lin1<=chr2[j,3]){

cat(chr1[i,],"\n",append = TRUE,file="snpsingene.txt")

}

    }

}

(#耗时较长，改进策略1'将命令拆分成多个任务同时进行；2'直接对基因组注释的位置进行操作比对；)