ABSOLUTE肿瘤纯度分析测试

library(data.table)
library(foreach)
library(doMC)
library(ABSOLUTE)

setwd("/home/train/data/ABSOLUTE_test/Tumor_purity")

DoAbsolute <- function(scan,input) {
#  registerDoSEQ()  #多核心计算时候,需要注释掉
  library(ABSOLUTE)
  plate.name <- "test" #可以根据需要修改
  genome <- "hg18"
  platform <- c("SNP_6.0", "Illumina_WES", "SNP_250K_STY")
  sample.name <- scan
  primary.disease = "NA"
  sigma.p <- 0
  max.sigma.h <- 0.02
  min.ploidy <- 0.95
  max.ploidy <- 10
  max.as.seg.count <- 1500
  max.non.clonal <- 0
  max.neg.genome <- 0
  copy_num_type <- "total"
  seg.dat.fn <- input
  results.dir <- file.path(".", "output", scan, "absolute")
  print(paste("Starting scan", scan, "at", results.dir))
  log.dir <- file.path(".", "output", "abs_logs")
  if (!file.exists(log.dir)) {
     dir.create(log.dir, recursive=TRUE)
  }
  dz <- lapply(results.dir,function(x){
    if(!file.exists(x)){
      dir.create(x, recursive=TRUE)
    }
  })
  rm(dz)
  sink(file=file.path(log.dir, paste(scan, ".abs.out.txt", sep="")))
  RunAbsolute(seg.dat.fn, sigma.p, max.sigma.h, min.ploidy, max.ploidy, primary.disease, 
              platform, sample.name, results.dir, max.as.seg.count, max.non.clonal, 
              max.neg.genome, copy_num_type, verbose=TRUE)
  sink()
}


registerDoMC(20)    #对每组数据调用线程池

arrays.txt <- "SNP6_solid_tumor.seg.txt"

input_data <- "test.txt"
scans <- unique(as.data.frame(fread(arrays.txt))$Sample)
foreach (scan=scans, .combine=c) %dopar% {
  DoAbsolute(scan,input=input_data)
}

obj.name <- "test_summary"
results.dir <- file.path(".", "output", "abs_summary")
absolute.files <- file.path(".", "output",
                            scans, "absolute",
                            paste(scans, ".ABSOLUTE.RData", sep=""))


CreateReviewObject(obj.name, absolute.files, results.dir, "total", verbose=TRUE)

calls.path = file.path("output", "abs_summary", "test_summary.PP-calls_tab.txt")
modes.path = file.path("output", "abs_summary", "test_summary.PP-modes.data.RData")
output.path = file.path("output", "abs_extract")
ExtractReviewedResults(calls.path, "test", modes.path, output.path, "absolute", "total")


# doabsolute包更好使用,需要根据输入数据调整格式暂时分步处理,今后有时间写成pipline可以无缝对接doabsolute进行批量处理

mkdir T1_5
mkdir total
mv T1_5.cr.igv.seg T1_5
mv T1_5_absolute.maf T1_5
cd T1_5/
mv T1_5.cr.igv.seg tumor.seg
mv ../T1_3/step1_Xto23.sh .
bash step1_Xto23.sh 
cp ../T1_3/step2_MafFilebalance.pl .
perl step2_MafFilebalance.pl T1_5_absolute.maf
cp ../T1_3/step3_SegFilebalance.pl .
perl step3_SegFilebalance.pl new_tumor.seg
cp ../T1_3/step4_rebalance_seg_and_maf.R .
Rscript step4_rebalance_seg_and_maf.R
cp ../T1_3/step5_get_seperate_file.R .
Rscript step5_get_seperate_file.R
cp ../T1_3/step6_prepare_doabsolute.pl .
perl step6_prepare_doabsolute.pl ./seperate_file/
cp ../T1_3/step6.5_prepare_doabsolute.sh .
#修改step6.5中的文件名为T1_5
bash step6.5_prepare_doabsolute.sh

cd ..
mkdir total
cd total/
mkdir seperate_file
cp ../T1_3/seperate_file/*.maf seperate_file/
cp ../T1_3/seperate_file/new_T1_3.seg seperate_file/
cp ../T1_4/seperate_file/*.maf seperate_file/
cp ../T1_4/seperate_file/new_T1_4.seg seperate_file/
cp ../T1_5/seperate_file/*.maf seperate_file/
cp ../T1_5/seperate_file/new_T1_5.seg seperate_file/
cp ../T1_3/new_final.seg .
cat ../T1_4/new_final.seg | awk 'NR>1' >> new_final.seg
cat ../T1_5/new_final.seg | awk 'NR>1' >> new_final.seg
cp ../T1_3/step7_doabsolute.R .
#注意修改路径
Rscript step7_doabsolute.R
发表在 R | 留下评论

R 函数构造练习

name <- rownames(rt_expr)   #基因名

immune <- function(genename){
    y <- as.numeric(rt_expr)    #单一基因的表达量表格
    do.call(rbind,lapply(colnames(rt),function(x){    #rt列明与rt_expr的行名相同 
        test <- cor.test(rt[,x],y,type="spearman")
        cor <- test$estimate
        p.value <- test$p.value
        data.frame(gene=genename,symbol=x,correlation=cor,pvalue=p.value)
        }))
    }

data <- immune(name)


函数的参数的传入也可以多个,以及固定参数如下: 
compute <- function(x,y,d=2){
    ....
    }
发表在 R | 一条评论

perl: warning: Falling back to the standard locale (“C”)

Win10子系统登录服务器,提示本地电脑缺少相应的语言库,并报错perl: warning: Falling back to the standard locale (“C”),解决方式如下:

sudo apt-get install language-pack-en
sudo locale-gen
sudo vim /etc/default/locale
#添加
    LC_ALL=en_US.utf8 
cd /usr/share/locales/
sudo ./install-language-pack en_US 
再重启终端即可
发表在 Linux | 留下评论

C++对二进制和文本文件操作练习

#文本文件

#include <iostream>
#include <fstream>

using namespace std;

int main()
{
    int a[10],max,i,order;
    ifstream infile("f1.dat",ios::in);
    {
        cout << "open error !" << endl;
        exit(1);
    }
    for (i=0;i<10;i++)
    {
        infile >> a[i];
        cout << a[i] << " ";
    }
    cout << endl;
    max=a[0];
    order=0;
    for (i=1;i<10;i++)
    {
        if(a[i]>max)
        {
            max=a[i];
            order=i;
        }
    }
    cout << "max=" << max << endl << "order=" << order << endl;
    infile.close();

    return 0;
} 

#二进制文件

#include <fstream>
#include <iostream>

using namespace std;

struct student
{
    int num;
    char name[20];
    float score;
};

int main()
{
    student stud[5]={1001,"Li",85,1002,"Fun",97.5,1004,"Wang",54,1006,"Tan",75.6,1010,"ling",96};
    fstream iofile("stud.dat",ios::in|ios::out|ios::binary);
    if(!iofile)
    {
        cerr << "open error !" << endl;
        abort();
    }
    for(int i=0;i<5;i++)
        iofile.write((char*)&stud[i],sizeof(stud[i]));
    student stud1[5];
    for (int i=0;i<5;i=i+2)
    {
        iofile.seekg(i*sizeof(stud1[i]),ios::beg);
        iofile.read((char*)&stud1[i/2],sizeof(stud1[i/2]));
        cout << stud1[i/2].num << " " << stud1[i/2].name << " " << stud1[i/2].score << endl;
    }
    cout << endl;
    stud[2].num=1012;
    strcpy(stud[2].name,"Wu");
    stud[2].score=60;
    iofile.seekp(2*sizeof(stud[0]),ios::beg);
    iofile.write((char*)&stud[2],sizeof(stud[2]));
    iofile.seekg(0,ios::beg);
    for(int i=0;i<5;i++)
    {
        iofile.read((char*)&stud[i],sizeof(stud[i]));
        cout << stud[i].num << " " << stud[i].name << " " << stud[i].score << endl;
    }
    iofile.close();

    return 0;
}
发表在 C++ | 留下评论

虚函数及捕获异常练习

#虚函数
#include <iostream>
using namespace std;

class base  //定义基类
{
    public:
        virtual void who()  //虚函数声明
        {
            cout<<"this is the class of base !"<<endl;
        }
};

class derive1: public base  //定义派生类derive1
{
    public:
        void who()  //重新定义虚函数
        {
            cout << "this is the class of derive1 !" << endl;
        }
};

class derive2: public base  //定义派生类derive2
{
    public:
        void who()  //重新定义虚函数
        {
            cout << "this is the class of derive2 !" << endl;
        }
};

int main()
{
    base obj, *ptr; //声明基类对象obj、指针ptr
    derive1 obj1;   //声明派生类1的对象obj1
    derive2 obj2;   //声明派生类2的对象obj2
    ptr= &obj;  //基类指针指向基类对象
    ptr -> who();   //调用基类成员函数
    ptr = &obj1;    //基类指针指向派生类1对象
    ptr -> who();   //调用派生类1成员函数
    ptr=&obj2;  //基类指针指向派生类2对象
    ptr -> who();   //调用派生类2成员函数
    return 0;
}

#捕获异常
#include <iostream>

using namespace std;

int Div(int x,int y);

int main()
{
    try
    {
        cout<<"5/2="<<Div(5,2)<<endl;
        cout<<"8/0="<<Div(8,0)<<endl;
        cout<<"7/1="<<Div(7,1)<<endl;   //由于8/0抛出异常,所以7/1并没有运行
    }
    catch(int)
    {
        cout<<"except of deviding zero"<< endl;;
    }
    cout << "that is ok" << endl;
    return 0;
}

int Div(int x,int y)
{
    if(y==0)
    {
        throw y;
    }
    return x/y;
}
发表在 C++ | 留下评论

基类对象与派生类对象指针

#include <iostream>
#include <string>

using namespace std;

class String
{
    char *name;
    int length;
    public:
        String(char *str)
        {
            length=strlen(str);
            name=new char[length+1];
            strcpy(name,str);
        }
        void show()
        {
            cout<<name<<endl;
        }
};

class de_string:public String
{
int age;
public:
    de_string(char *str, int age):String(str)
    {
        de_string::age=age;
    }
    void show()
    {
        String::show();
        cout<<"the age is:"<<age<<endl;
    }
};

int main()
{
    String s1("tony"),*ptr1;
    de_string s2("jack",20),*ptr2;
    ptr1=&s1;
    ptr1->show();
    ptr1=&s2;   //将ptr1指向String类的派生类de_string的对象s2
    ptr1->show();   //调用s2对象所属的基类的成员函数show()
    ptr2=&s2;
    ptr2 -> show();
    return 0;
}
发表在 C++ | 留下评论