scDNA-seq genomic analysis pipline

 

数据下载

这次用到的单细胞数据来自于文献 Multiclonal Invasion in Breast Tumors Identified by Topographic Single Cell Sequencing. 作者并没有公布处理完的数据,只把原始测序数据上传到了 NCBI 的 SRA 数据库, 一共 10 个 synchronous ductal carcinoma.
为了方便处理,我们将数据按照样本来源分别下载存储,在下载时遵循以下顺序:
1. 确认样本来源以及基础信息

  1. 进入样本主页,浏览数据信息

  2. 确认每个细胞数据来源,下载 SraRunInfo 文件

由于样本数据过多, wget 和 axel　在下载速度上都有所缺陷，需要用迅雷批量下载，这就需要首先产生所有需要下载样本的 URL　地址: 例如样本 SRR6238488 存储的 URL 为 : ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR623/SRR6238488/SRR6238488.sra

setwd("H://P1")
cell_list <- unlist(read.csv("SraRunInfo.csv", header = T, as.is = T)[,1])
header <- "ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/"
url <- matrix(apply(
    matrix(cell_list, ncol = 1),
    1,
    function(x){
        folder = paste(unlist(strsplit(x, ""))[1:6], collapse = "")
        url = paste(header, folder, "/", x, "/", x, ".sra", sep = "")
        return(url)
    }
), ncol = 1)
write.table(url,"url.txt", col.name = F, row.names = F, quote = F)

下载完一个样本的所有细胞　sra　数据后，就可以用 fastq-dump处理得到对应的 fasta 文件了.

cell_list <- (read.csv(
    "/media/wang/WWD/P1/SraRunInfo.csv",
    header = T,
    stringsAsFactors = F
    )[,1]) %>%  unlist()
folder <- "/media/wang/WWD/P1/cell"
for(cell in cell_list){
  str = paste0(
    "fastq-dump -O /media/wang/WWD/P6/cell/fasta",
    "--split-3 /media/wang/WWD/P6/cell/", cell, ".sra"
    )
  system(str)
}

tophat进行序列比对

利用 Bowtie2 构建参考序列 index:

bowtie2-build
/media/wang/subclone_anlysis/data/ref/hg19.fa
/media/wang/subclone_anlysis/data/ref/hg19

一开始是打算用 STAR 进行参考序列构建和序列比对的，但是它瞧不起我，我不开心了，我有脾气了.

接下去就是批量跑 TopHat:

cell_list <- (read.csv(
    "/media/wang/WWD/P1/SraRunInfo.csv",
    header = T,
    stringsAsFactors = F
    )[,1]) %>%  unlist()

for(cell in cell_list)){
  str = paste0(
    "tophat -o /media/wang/Seagate/subclone/bam",
    "-G /media/wang/subclone_anlysis/data/ref/GRCh37.gtf",
    "/media/wang/subclone_anlysis/data/ref/hg19",
    "/media/wang/subclone_anlysis/data/sra/cell/p9/fasta/", cell, ".fastq"
    )
  system(str)
  str = paste0(
    "samtools rmdup -s /media/wang/subclone/bam/accepted_hits.bam",
    "/media/wang/Seagate/subclone/bam/", cell, ".bam"
    )
  system(str)
}

进行基因组分析

得到了每个细胞对应的 bam 文件以后，利用 varscan工具进行 germline variants calling 以及 Copy Number Alteration (CNA) Calling.
由于前面的tophat结果还没跑完，所以只有两个示例:

  1. germline variants calling
samtools mpileup -q 1 -f <BR>
/media/wang/subclone_anlysis/data/ref/hg19.fa normal.bam tumor.bam | <BR>
java -jar VarScan.v2.3.9.jar somatic --output-file CNA --mpileup 1
  1. CNA calling
samtools mpileup -q 1 -f
/media/wang/subclone_anlysis/data/ref/hg19.fa normal.bam tumor.bam |
java -jar VarScan.v2.3.9.jar copynumber --output-file CNA --mpileup 1

以上就是关于单细胞测序数据的完整基因组分析数据处理流程，完结撒花.

posted @ 2018-12-13 20:00  PeRl`  阅读(585)  评论(2)  编辑  收藏  举报