r3fang / SnapATAC

Analysis Pipeline for Single Cell ATAC-seq
GNU General Public License v3.0
296 stars 124 forks source link

"error: the following nsap files do not contain AM session" #177

Closed rojinsafavi closed 4 years ago

rojinsafavi commented 4 years ago

I have five atac files. I ran snaptools for each one of them: snaptools snap-pre --input-file=fragments.bed.gz --output-snap=fragments.snap --genome-name=hg38 --genome-size=../../hg38.chrom.sizes

then snaptools snap-add-bmat for each one of them: snaptools snap-add-bmat --snap-file=fragments.snap --bin-size-list 1000 5000 10000 --verbose=True

then:

library(SnapATAC);
snap.files = c(
    "one/fragments.snap", 
    "two/fragments.snap",
    "three/fragments.snap",
    "four/fragments.snap",
    "five/fragments.snap"
  );
sample.names = c(
    "one",
    "two",
    "three",
    "four",
    "five"
  );
barcode.files = c(
    "one/singlecell.csv",
    "two/singlecell.csv",
    "three/singlecell.csv",
    "four/singlecell.csv",
    "five/singlecell.csv"
  );
x.sp.ls = lapply(seq(snap.files), function(i){
    createSnap(
        file=snap.files[i],
        sample=sample.names[i]
    );
  })
names(x.sp.ls) = sample.names;
barcode.ls = lapply(seq(snap.files), function(i){
    barcodes = read.csv(
        barcode.files[i], 
        head=TRUE
    );
    # remove NO BAROCDE line
    barcodes = barcodes[2:nrow(barcodes),];
    barcodes$logUMI = log10(barcodes$passed_filters + 1);
    barcodes$promoter_ratio = (barcodes$promoter_region_fragments+1) / (barcodes$passed_filters + 1);
    barcodes
  })

Then:

x.sp.ls

number of barcodes: 461144
number of bins: 0
number of genes: 0
number of peaks: 0
number of motifs: 0

$two
number of barcodes: 468172
number of bins: 0
number of genes: 0
number of peaks: 0
number of motifs: 0

$three
number of barcodes: 443411
number of bins: 0
number of genes: 0
number of peaks: 0
number of motifs: 0

$four
number of barcodes: 476659
number of bins: 0
number of genes: 0
number of peaks: 0
number of motifs: 0

$five
number of barcodes: 471234
number of bins: 0
number of genes: 0
number of peaks: 0
number of motifs: 0```

Then:

for both datasets, we identify usable barcodes using [3.5-5] for log10(UMI) and [0.4-0.8] for promoter ratio as cutoff.

cutoff.logUMI.low = c(3.5, 3.5, 3.5, 3.5, 3.5); cutoff.logUMI.high = c(5, 5, 5, 5, 5); cutoff.FRIP.low = c(0.4, 0.4, 0.4, 0.4, 0.4); cutoff.FRIP.high = c(0.8, 0.8, 0.8, 0.8, 0.8);

barcode.ls = lapply(seq(snap.files), function(i){ barcodes = barcode.ls[[i]]; idx = which( barcodes$logUMI >= cutoff.logUMI.low[i] & barcodes$logUMI <= cutoff.logUMI.high[i] & barcodes$promoter_ratio >= cutoff.FRIP.low[i] & barcodes$promoter_ratio <= cutoff.FRIP.high[i] ); barcodes[idx,] });

x.sp.ls = lapply(seq(snap.files), function(i){ barcodes = barcode.ls[[i]]; x.sp = x.sp.ls[[i]]; barcode.shared = intersect(x.sp@barcode, barcodes$barcode); x.sp = x.sp[match(barcode.shared, x.sp@barcode),]; barcodes = barcodes[match(barcode.shared, barcodes$barcode),]; x.sp@metaData = barcodes; x.sp }) names(x.sp.ls) = sample.names; x.sp.ls

$one number of barcodes: 1794 number of bins: 0 number of genes: 0 number of peaks: 0 number of motifs: 0

$two number of barcodes: 1077 number of bins: 0 number of genes: 0 number of peaks: 0 number of motifs: 0

$three number of barcodes: 1886 number of bins: 0 number of genes: 0 number of peaks: 0 number of motifs: 0

$four number of barcodes: 1077 number of bins: 0 number of genes: 0 number of peaks: 0 number of motifs: 0

$five number of barcodes: 2518 number of bins: 0 number of genes: 0 number of peaks: 0 number of motifs: 0

x.sp = Reduce(snapRbind, x.sp.ls); x.sp@metaData["sample"] = x.sp@sample; x.sp

number of barcodes: 8352 number of bins: 0 number of genes: 0 number of peaks: 0 number of motifs: 0


```table(x.sp@sample);```

onw two three four five 1794 2518 1886 1077 1077


And I get this error here:

x.sp = addBmatToSnap(x.sp, bin.size=5000)```

[1] "error: the following nsap files do not contain AM session"
[[1]]
[1] "three/fragments.snap"

[[2]]
[1] "five/fragments.snap"

Error in addBmatToSnap.default(x.sp, bin.size = 5000): 
Traceback:

1. addBmatToSnap(x.sp, bin.size = 5000)
2. addBmatToSnap.default(x.sp, bin.size = 5000)
3. stop()