数据下载
这次用到的单细胞数据来自于文献 Multiclonal Invasion in Breast Tumors Identified by Topographic Single Cell Sequencing. 作者并没有公布处理完的数据,只把原始测序数据上传到了 NCBI 的 SRA 数据库, 一共 10 个 synchronous ductal carcinoma. 
为了方便处理,我们将数据按照样本来源分别下载存储,在下载时遵循以下顺序:
1. 确认样本来源以及基础信息 
进入样本主页,浏览数据信息
确认每个细胞数据来源,下载
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结果还没跑完，所以只有两个示例:
- 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- 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以上就是关于单细胞测序数据的完整基因组分析数据处理流程，完结撒花.
