SNS: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
Line 24: Line 24:
*This method was developed by Peter Park's group at Harvard. [[http://www.pnas.org/content/108/46/E1128.short Xi et al. PNAS 2011]]
*This method was developed by Peter Park's group at Harvard. [[http://www.pnas.org/content/108/46/E1128.short Xi et al. PNAS 2011]]
*It starts with mapped bam files after removal of clonal reads.
*It starts with mapped bam files after removal of clonal reads.
*Step 1: generate read count files from the bam file (replace DIR with a directory name to write all the .seq files).
*Step 1: generate read count files from the bam file (<b>replace DIR with a directory name to write all the .seq files</b>).
   /home/kunzhang/softwares/BICseq/PERL/BICseq_1.1.2/SAMgetUnique/samtools-0.1.7a_getUnique-0.1.1 view -U BWA,DIR/,N,N SAMPLE_NAME.bam
   /home/kunzhang/softwares/BICseq/PERL/BICseq_1.1.2/SAMgetUnique/samtools-0.1.7a_getUnique-0.1.1 view -U BWA,DIR/,N,N SAMPLE_NAME.bam
   Repeat the same procedure for each bam files, use one output directory per sample
   Repeat the same procedure for each bam files, use one output directory per sample

Revision as of 06:23, 15 September 2012

CNV calling on single-cell genome sequencing data

CSHL SNS caller

  • The method was developed by CSHL Baslan et al. Nat Protoc, 2012.
  • A Perl script fastq2cnv.pl was written to implement this method in both Triton and Genome-miner.
  • To run this script, one needs a .info file per sample. The file should contain the server name (the first line, either triton or genome-miner), location and other information of the sequencing data. Here is an example of a .info file.
  • On genome-miner, the job can be started with the following command:
  nohup PATH_NAME/fastq2cnv.pl  SAMPLE_NAME.info > SAMPLE_NAME_SNS.log &
  • On Triton, a job file needs to be created for each sample. And the job can be submitted by qsub. Here is an example of .job file.
  qsub SAMPLE_NAME.job
  • Once the job is completed, transfer the SAMPLE_NAME.varbin50k.out to your local computer (you can also run R in genome-miner, but not triton).
  • Edit two R scripts by replacing the sample name with the actual one you are using.
  SAMPLE_NAME.cbs.r
  SAMPLE_NAME.copynumber.r
  • Under R, run the following two commands (make sure both the varbin50k.out file and the R scripts are in the current directory).
  source("SAMPLE_NAME.cbs.r");
  source("SAMPLE_NAME.copynumber.r");
  • You will find two dozens of plots in the Postscript format, plus a number of other text files. Postscript files can be converted into pdf files using Acrobat Distiller, or ps2pdf in Unix/Linux.
   SAMPLE_NAME.hg19.50k.k50.varbin.data.copynumber.
   SAMPLE_NAME.hg19.50k.k50.nobad.varbin.data
   SAMPLE_NAME.hg19.50k.k50.nobad.varbin.short
   SAMPLE_NAME.hg19.50k.k50.varbin.data
   SAMPLE_NAME.hg19.50k.k50.varbin.short

BIC-Seq

  • This method was developed by Peter Park's group at Harvard. [Xi et al. PNAS 2011]
  • It starts with mapped bam files after removal of clonal reads.
  • Step 1: generate read count files from the bam file (replace DIR with a directory name to write all the .seq files).
 /home/kunzhang/softwares/BICseq/PERL/BICseq_1.1.2/SAMgetUnique/samtools-0.1.7a_getUnique-0.1.1 view -U BWA,DIR/,N,N SAMPLE_NAME.bam
 Repeat the same procedure for each bam files, use one output directory per sample
  • Step 2: create a tab-delimited config file as the following:
 chrom   case    control
 chr1    case/chr1.seq   control/chr1.seq
 chr2    case/chr2.seq   control/chr2.seq
  • Step 3: run BIC-seq.pl:
 /home/kunzhang/softwares/BICseq/PERL/BICseq_1.1.2/BIC-seq/BIC-seq.pl [options] <ConfigFile> <OutputDir> <Description>
 Options:
       --help
       --lambda=<float>: default 2
       --bin_size=<int>: default 100
       --multiplicity=<float>: default 2
       --window=<int>: the window for removing the outliers; default 200
       --f=<float>: expected number of type I errors in the merging process; An alternative way to specify lambda
       --B=<int>: number of permutations for FDR estimate; default 0
       --paired: if specified the data is treated as paired-end data
       --I=<Insert,SDofInsert>: specify the insert size and standard deviation of insert size. Default <200,20>
 <Outputdir>: the name of the directory to write the output (this cannot be an existing directory).
 <Description>: used to name the segmentation results under <OutputDir>
  • For detail, see /home/kunzhang/softwares/BICseq/PERL/BICseq_1.1.2/BIC-seq/README