I'm trying to use Seurat for the first time and am learning about single-cell analysis for the first time, and I'm doing so with split-seq data. (Full disclosure: I'm also a lightweight when it comes to R.)
When Read10X didn't work on the mtx and associated files output by the split-seq-pipeline no matter how much I tried to make it look exactly like a 10x dataset (based on numerous examples and descriptions), I tried converting my sparse market matrix format into a dense tsv format (see update below for a sample of the file) and I read it in using:
dge_data <- as.matrix(read.table("DGE.tsv", sep="\t", header=TRUE)) pbmc <- CreateSeuratObject(counts = dge_data, min.cells = 3, min.features = 200, project = "SPLITSEQ_TT", assay = "RNA")
I started following chapter 9 of the scRNA-Seq course I found here. I ran into some issues which I was able to resolve by tweaking my mtx to tsv conversion script, such as including the chromosome in the row names (so that the mitochondrial QC step would work), but I ran into some errors that I assume have to do with the format of the gene names/IDs and possibly some issues stemming from cell names/IDs (which I've thus far been able to work through I think).
However, given errors like the following, and the fact that I can't seem to find specifications on these sorts of row/column name requirements (e.g. no commas and must include at least the "MT" chromosome name for QC steps using mitochondrial genes), I was wondering if I could be directed to somewhere where those specifications are laid bare. For the most part, the gene names have not been a problem, but I run afoul in this step (where it appears to complain about a comma):
> FeatureScatter(object = pbmc, feature1 = "hg38-4,ENSG00000074966,TXK", feature2 = "PC_1") Error in parse(text = x) : <text>:1:7: unexpected ',' 1: hg38.4, ^
Also, when I look at how CellsByIdentities behaves (while I was working through other errors/warnings), I realized I'm not sure that the column names are being processed correctly due to the underscores (followed by a number (e.g.
_0)), because my understanding is that the
_# portion of the split-seq cell barcode indicates the first well barcode, which would be needed to uniquely identify the cells in the table. I had retrieved the cell barcode IDs from the cell_metadata.csv file output by split-seq. Should I modify these column names? Because I infer from this output from CellsByIdentities that it's only using the sequence portion of the cell barcodes to uniquely identify cells, which I don't think is correct:
> CellsByIdentities(object=pbmc, cells=colnames(pbmc)) $AAACATCGAAACATCG  "AAACATCGAAACATCG_0" "AAACATCGAAACATCG_1" $AAACATCGAACGTGAT  "AAACATCGAACGTGAT_0" $AACGTGATAAACATCG  "AACGTGATAAACATCG_0" "AACGTGATAAACATCG_1" $AACGTGATAACGTGAT  "AACGTGATAACGTGAT_0" "AACGTGATAACGTGAT_1"
The other assumptions have to do with data size. I was running a small set of test data, and I had to adjust some parameters. So I would like to know if there's a hard lower limit on the data. And do I need to restructure my data to separate counts of genes from different species (e.g. mouse & human)? I.e. Should I create separate matrices that have only the human-identified cells and gene annotations to use it with Seurat?
I can't seem to find these specifications, presumably because Seurat was designed exclusively around 10x data and that the tool wasn't tested on data from other sources. Perhaps I'm worrying too much, but without knowing what may rely on column/row name format assumptions, I can't be sure my data is being analyzed properly.
Here is the head of the tsv file I originally tried to load when I ran into various errors that I suspect were due to the formatting of the row and column names (note, there are only 7 cells total, and it's a small but real test dataset from our core facility):
AAACATCGAAACATCG_0 AAACATCGAAACATCG_1 AAACATCGAACGTGAT_0 AACGTGATAAACATCG_0 AACGTGATAAACATCG_1 AACGTGATAACGTGAT_0 AACGTGATAACGTGAT_1 hg38_X,ENSG00000000003,TSPAN6 8 2 2 26 15 17 13 hg38_X,ENSG00000000005,TNMD 0 0 0 0 0 2 0 hg38_20,ENSG00000000419,DPM1 22 26 27 73 43 50 28 hg38_1,ENSG00000000457,SCYL3 6 6 5 12 7 19 8 hg38_1,ENSG00000000460,C1orf112 23 10 24 59 34 49 21 hg38_1,ENSG00000000971,CFH 101 44 129 331 151 270 122 hg38_6,ENSG00000001036,FUCA2 33 28 24 82 58 74 36 hg38_6,ENSG00000001084,GCLC 30 16 14 72 27 54 42 hg38_6,ENSG00000001167,NFYA 14 9 17 41 28 31 12
Yesterday, I did end up trying to change the column and row names to the following:
AAACATCGAAACATCG-0 AAACATCGAAACATCG-1 AAACATCGAACGTGAT-0 AACGTGATAAACATCG-0 AACGTGATAAACATCG-1 AACGTGATAACGTGAT-0 AACGTGATAACGTGAT-1 hg38-X.ENSG00000000003.TSPAN6 8 2 2 26 15 17 13 hg38-X.ENSG00000000005.TNMD 0 0 0 0 0 2 0 hg38-20.ENSG00000000419.DPM1 22 26 27 73 43 50 28 hg38-1.ENSG00000000457.SCYL3 6 6 5 12 7 19 8 hg38-1.ENSG00000000460.C1orf112 23 10 24 59 34 49 21 hg38-1.ENSG00000000971.CFH 101 44 129 331 151 270 122 hg38-6.ENSG00000001036.FUCA2 33 28 24 82 58 74 36 hg38-6.ENSG00000001084.GCLC 30 16 14 72 27 54 42 hg38-6.ENSG00000001167.NFYA 14 9 17 41 28 31 12
And it seems to have fixed most of the issues. The violin plots now have curvy brown swaths down the center that were not there before. And in particular, when I run CellsByIdentities manually on the column names now, I get what looks correct (as opposed to what I pasted yesterday above (which appeared to group the 7 cells in 3 pairs and one individual):
CellsByIdentities(object = pbmc, cells = colnames(pbmc)) $SPLITSEQ_TT  "AAACATCGAAACATCG.0" "AAACATCGAAACATCG.1" "AAACATCGAACGTGAT.0" "AACGTGATAAACATCG.0" "AACGTGATAAACATCG.1" "AACGTGATAACGTGAT.0"  "AACGTGATAACGTGAT.1"
I think that the remaining issues are likely due to the small data size. For reference as to the data size, this is the first data line from the original mtx file:
7 38526 194539
Please note also that I'm not looking for how to get this particular dataset to work. I am (as a member of the bioinformatics team that supports users of our core sequencing facility) trying to establish the requirements for the conversion script so that when we deliver SPLiT-Seq results, our output can be directly used with Seurat by the users without issues or need for further conversion. So I specifically need to know what makes row or column names valid or invalid, as opposed to what to change WRT this specific dataset. We may be supplied with annotations from various sources that may have different characters in the gene symbols, or accession numbers, or even chromosome names.