text:
Cell Ranger) Error during performing a test run
I am interested in the analysis of single-cell RNA-sequencing data.
I installed Cell Ranger according to the instruction (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/tutorial_in#download).
I know, it should be easy, but I found an error occurring after test running as below.
[error] Pipestance failed. Error log at:
tiny/SC_RNA_COUNTER_CS/SC_RNA_COUNTER/_BASIC_SC_RNA_COUNTER/COUNT_GENES/fork0/chnk0-u6c389ca928/_errors

Within that file, it shows the following error:
Job failed in stage code. signal: killed.

So, what do you think about this problem?
Additional information:

a checksum test (to check the integrity of the downloaded files) was successful
I performed a test run on a personal laptop, mac.
I used SSH command to connect the server, but we are sharing the server with our lab members.

I found an error as follows:
WARNING: configured to use 112GB of local memory, but only 80.5GB is currently available.

However, at this time with unknown reasons, the installation is successfully completed, in spite of the warning. Do you have idea?
Thank you in advance.

I think I am having the same issue. but the error occurred while running cellranger-arc. If I am sharing server (hpc) with other people how do i solve this problem ? Thanks

Hello @Goku, if you have a problem that is similar to this, but not really, open a new question, link this one and explain properly why this solution did not work for you. Posting answers as follow up questions is not the way of stack exchange :-).

Answer from PPK, converted from comments:
The memory warning is only letting you know that the upper limit you set on memory (112 gb) was higher than the available memory. Do you share the server with other people? because the minimum amount of memory for Cell Ranger is 64 GB. In your initial run there may have been less than 64 GB free.


text:
What does it mean when PDB lists multiple organisms for the same structure?
I'm browsing through ribosomes which have been modeled with Cryoem on PDB and am pretty confused by the fact that for some, PDB, or rather the authors of the deposition probably, put multiple organisms from different domains of life... How does one interpret that?
So far as i understand bos taurus is a cow and e.coli is a bacterium. So eukarya and bacteria.. Thank you.

Since this is a 40S ribosomal protein, and therefore eukaryotic-specific (as the description says) the listing of E. coli is not intuitive as you note. If you go down to "Macromolecules" and click on "Nucleic Acids/Hybrids" you can see the listed E. coli sequence is just 28bp, and all Us (poly(U) tail?). So, you can safely ignore the bacterial annotation I think.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
How to concatename file names in Snakemake?
I want to pass a series of files (e.g., a, b) to a program prog like: prog a,b.
Is there a way to do this in Snakemake?
My attempt is:
rule run_prog:
  input: "a", "b"
  output: "output"
  run: "prog {input}"

But the actual command Snakemake generated is prog a b instead of prog a,b.
Any solution to pass a comma concatenated file names to Snakemake?

run: "prog {input[0]},{input[1]}" ?

@Chris_Rands that works, is there a general solution, as there can be ~20 input files. Thanks!

You can use join to convert the list of input file into a string.
I think this would work as expected.:
rule run_prog:
    input: "a", "b"
    output: "output"
    run: f"prog {','.join(input)}"

This works great. It works for f"prog {','.join(input.variable)}"  if I define variable under input: section. Can you explain why have to use input[:] instead of input?

Right, using input by itself is the same ! I edited my answer.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
No MQ tags in VCF files
To call minority variants in my Mtb sequences I'm using a pipeline of
fastq ->
    BWA MEM for mapping ->
    picard MarkDuplicates ->
    GATK3 for indel realignment and base recalibration ->
    lofreq for variant calling (https://csb5.github.io/lofreq/)

However the resulting VCF files don't have an MQ tag for me to filter on. Where in this process would mapping quality information be generated or lost so that it doesn't end up in the VCF file?

MQ is a metric that characterises the alignment quality of an individual read. VCF contains variant calls for a base/region that is the product of all intersecting reads that contribute to the call being either REF or ALT. Therefore it makes no sense to report MQ in VCF. MQ is something that should be filtered for during variant calling, so the variant caller should have an option to only include alignments > MQ.

In GATK output for vcf, i think there are metrics that make use of MQ, for example MappingQuality, MappingQualityZero, MappingQualityRankSumTest. you can check it here https://gatk.broadinstitute.org/hc/en-us/articles/360037425531--Tool-Documentation-Index#VariantAnnotations

As @ATpoint explains in his comment, the mapping quality (MQ) applies to read mappings and, hence, is contained in BAM files, but not VCFs. The MQ values carry important information that ideally the variant caller should use when determining the variant quality (QUAL in the VCF) of individual variant loci.
In the LoFreq documentation it says that

Choice of mapper: Bowtie and BWA-SW
BWA-SW assigns very low mapping qualities to mapped reads. LoFreq incorporates mapping quality by default into its variant calling model, which doesn’t help in this scenario and will reduce its sensitivity. Consider disabling the use of mapping quality for lofreq call with -N.

So LoFreq is already aware of the MQ values and takes them into consideration during variant calling.
Since you are using BWA-MEM (and not BWA-SW) you should probably not be concerned about the issue they mention and trust their quality filtering.
If you are sceptic about the quality of mappings, though, you could still filter them out directly after mapping, e.g. using samtools view -q 20 ...


text:
Adding certain barcodes to metadata in Seurat
I want to upload an excel file sheet that has certain barcodes that I would like to show on my umap. How do I go about adding the file and linking it to the metadata?
Below is my following code.
# set up the working directory
wd = "/home/PTX_AAC656
setwd(wd)

# load counts
# SD-52-JZ33 patient sample but aligned to both mouse and human
SD_52_JZ33.data <- Read10X(data.dir = "../PTX_AAC656/metadata/patient_sample/SD-52-JZ33/outs/filtered_feature_bc_matrix/")
BNM_human_only.data <- Read10X(data.dir = "../PTX_AAC656/metadata/human_only_alignment/SD-52-JZ33_BNM_5K/outs/filtered_feature_bc_matrix/")
BNL_human_only.data <- Read10X(data.dir = "../PTX_AAC656/metadata/human_only_alignment/SD-52-JZ33_BNL_5K/outs/filtered_feature_bc_matrix/")
BNL_mouse_only.data <- Read10X(data.dir = "../PTX_AAC656/metadata/mouse_only_alignment/SD-52-JZ33_BNL_5K/outs/filtered_feature_bc_matrix/")
BNM_mouse_only.data <- Read10X(data.dir = "../PTX_AAC656/metadata/mouse_only_alignment/SD-52-JZ33_BNM_5K/outs/filtered_feature_bc_matrix/")

   BNL_both.data <- Read10X(data.dir = "../AAC656/SD-52-JZ33_BNL_5K/outs/filtered_feature_bc_matrix/")
    BNM_both.data <- Read10X(data.dir = "../AAC656/SD-52-JZ33_BNM_5K/outs/filtered_feature_bc_matrix/")

setwd(wd)

# create the Seurat objects
SD_52_JZ33 <- CreateSeuratObject(counts = SD_52_JZ33.data, min.cells = 0, project = "human")
BNL_human_only <- CreateSeuratObject(counts = BNL_human_only.data, min.cells = 0, project = "human_alignment_ptx")
BNM_human_only <- CreateSeuratObject(counts = BNM_human_only.data, min.cells = 0, project = "human_alignment_ptx")
BNM_both_ptx <- CreateSeuratObject(counts = BNM_both.data, min.cells = 0, project = "human_mouse_alignment_ptx")
BNL_both_ptx <- CreateSeuratObject(counts = BNL_both.data, min.cells = 0, project = "human_mouse_alignment_ptx")
BNL_both_ptx$model = "mouse_human_alignment_ptx"
BNM_both_ptx$model = "mouse_human_alignment_ptx"
SD_52_JZ33$model = "patient_aligned_mouse_human"
BNM_human_only$model = "human_alignment_only_ptx"
BNL_human_only$model = "human_alignment_only_ptx"
BNL_both_ptx$model = "mouse_human_alignment_ptx"
BNM_both_ptx$model = "mouse_human_alignment_ptx"
BNM_mouse_only <- CreateSeuratObject(counts = BNM_mouse_only.data, min.cells = 0, project = "mouse_alignment_ptx")
BNL_mouse_only <- CreateSeuratObject(counts = BNL_mouse_only.data, min.cells = 0, project = "mouse_alignment_ptx")
BNL_mouse_only$model = "mouse_alignment_only_ptx"
BNM_mouse_only$model = "mouse_alignment_only_ptx"

# merge samples
mouse_human_patient <- merge(BNL_mouse_only, y = c(BNM_mouse_only, SD_52_JZ33), add.cell.ids = c("A", "B", "C"), project = "HTB2876")
human_ptx_human_patient <- merge(BNL_human_only, y = c(BNM_human_only, SD_52_JZ33), add.cell.ids = c("D", "E", "F"), project = "HTB2876")
ptx_human_patient.1 <- merge(BNL_both_ptx, y = c(BNM_both_ptx, SD_52_JZ33), add.cell.ids = c("G", "H", "I"), project = "HTB2876")

# mouse ptx model with patient
# mark the mito genes
mito.genes <- grep(pattern = "^MT-", x = rownames(x = mouse_human_patient), value = TRUE)
mouse_human_patient[["percent.mt"]] <- PercentageFeatureSet(mouse_human_patient, pattern = "^MT-")
mouse_human_patient[["log_nCount_RNA"]] <- log2(mouse_human_patient[["nCount_RNA"]]+1)

# remove cells with <200 RNA molecules, or >6000 molecules, or >30% mito
mouse_human_patient <- subset(mouse_human_patient, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 30)
mouse_human_patient <- NormalizeData(mouse_human_patient, normalization.method = "LogNormalize", scale.factor = 10000)
mouse_human_patient <- FindVariableFeatures(mouse_human_patient, selection.method = "vst", nfeatures = 2000)
all.genes.mouse <- rownames(mouse_human_patient)
mouse_human_patient <- ScaleData(mouse_human_patient, features = VariableFeatures(object = mouse_human_patient), vars.to.regress = c("nCount_RNA"))
mouse_human_patient <- RunPCA(mouse_human_patient, features = VariableFeatures(object = mouse_human_patient))
mouse_human_patient <- FindNeighbors(mouse_human_patient, dims = 1:20)
mouse_human_patient <- FindClusters(mouse_human_patient, resolution = 0.4)
mouse_human_patient <- RunUMAP(mouse_human_patient, dims = 1:20)
DimPlot(mouse_human_patient, reduction = "umap", group.by = "model")

# human only alignment pix with patient
# mark the mito genes
mito.genes <- grep(pattern = "^MT-", x = rownames(x = human_ptx_human_patient), value = TRUE)
human_ptx_human_patient[["percent.mt"]] <- PercentageFeatureSet(human_ptx_human_patient, pattern = "^MT-")
human_ptx_human_patient[["log_nCount_RNA"]] <- log2(human_ptx_human_patient[["nCount_RNA"]]+1)
human_ptx_human_patient <- subset(human_ptx_human_patient, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 30)

# remove cells with <200 RNA molecules, or >6000 molecules, or >30% mito
human_ptx_human_patient <- NormalizeData(human_ptx_human_patient, normalization.method = "LogNormalize", scale.factor = 10000)
human_ptx_human_patient <- FindVariableFeatures(human_ptx_human_patient, selection.method = "vst", nfeatures = 2000)
all.genes.human.ptx <- rownames(human_ptx_human_patient)
human_ptx_human_patient <- ScaleData(human_ptx_human_patient, features = VariableFeatures(object = human_ptx_human_patient), vars.to.regress = c("nCount_RNA"))
human_ptx_human_patient <- RunPCA(human_ptx_human_patient, features = VariableFeatures(object = human_ptx_human_patient))
DimPlot(human_ptx_human_patient, reduction = "pca")
human_ptx_human_patient <- FindNeighbors(human_ptx_human_patient, dims = 1:20)
human_ptx_human_patient <- FindClusters(human_ptx_human_patient, resolution = 0.4)
human_ptx_human_patient <- RunUMAP(human_ptx_human_patient, dims = 1:20)
DimPlot(human_ptx_human_patient, reduction = "umap", group.by = "model")

# ptx and human
# mark the mito genes
mito.genes.ptx_patient <- grep(pattern = "^MT-", x = rownames(x = ptx_human_patient.1), value = TRUE)
ptx_human_patient.1[["percent.mt"]] <- PercentageFeatureSet(ptx_human_patient.1, pattern = "^MT-")
ptx_human_patient.1[["log_nCount_RNA"]] <- log2(ptx_human_patient.1[["nCount_RNA"]]+1)

# remove cells with <200 RNA molecules, or >6000 molecules, or >30% mito
ptx_human_patient.1 <- subset(ptx_human_patient.1, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 30)
ptx_human_patient.1 <- NormalizeData(ptx_human_patient.1, normalization.method = "LogNormalize", scale.factor = 10000)
ptx_human_patient.1 <- FindVariableFeatures(ptx_human_patient.1, selection.method = "vst", nfeatures = 2000)
all.genes.both.patient.1 <- rownames(ptx_human_patient.1)
ptx_human_patient.1 <- ScaleData(ptx_human_patient.1, features = VariableFeatures(object = ptx_human_patient.1), vars.to.regress = c("nCount_RNA"))
ptx_human_patient.1 <- RunPCA(ptx_human_patient.1, features = VariableFeatures(object = ptx_human_patient.1))
ptx_human_patient.1 <- FindNeighbors(ptx_human_patient.1, dims = 1:20)
ptx_human_patient.1 <- FindClusters(ptx_human_patient.1, resolution = 0.4)
ptx_human_patient.1 <- RunUMAP(ptx_human_patient.1, dims = 1:20)

#run plots
DimPlot(ptx_human_patient.1, reduction = "umap", group.by = "model")
DimPlot(ptx_human_patient.1, reduction = "umap", label = "TRUE")
saveRDS(ptx_human_patient.1, file = "derived_data/PTX_human_2876_res04_6000g.rds")
DimPlot(ptx_human_patient.1, reduction = "umap", group.by = "model")
DimPlot(ptx_human_patient.1, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
ptx_human_patient.1$clustering.idents <- Idents(object = ptx_human_patient.1)
Idents(ptx_human_patient.1) <- "clustering.idents"
new.cluster.ids <- c("Human PTX Tumor", "PTX Fibroblast", "Human Fibroblast", "PTX Macrophage", "Human PTX Tumor", "Human Mitochondrial", "Human Macrophage", "Human Endothelial", "Human T-cells", "Human PTX Tumor","Human Fibroblast", "PTX Fibroblast", "Human Tumor Ductal", "Human Acinar", "PTX Macrophage", "PTX Stellate", "PTX Endothelial")
names(new.cluster.ids) <- levels(ptx_human_patient.1)
ptx_human_patient.1 <- RenameIdents(ptx_human_patient.1, new.cluster.ids)
saveRDS(ptx_human_patient.1, file = "derived_data/PTX_human_2876_res04_6000g.rds")
p1 = DimPlot(ptx_human_patient.1, reduction = "umap", label = TRUE, pt.size = 0.5,  group.by = "clustering.idents") + NoLegend()
p2 = DimPlot(ptx_human_patient.1, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
print(CombinePlots(plots = list(p1, p2)))
markers2 <- FindAllMarkers(ptx_human_patient.1, min.pct = 0.25, logfc.threshold = 0.25)
markers2 %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
saveRDS(markers, file="derived_data/markers_HTB2876_PC20_res04_6000g.rds")

Part of my excel file
barcode muTaTionsTaTus
GAGATGGCACTAGGCC-1  G12V
CTGAGCGAGACTCAAA-1  G12V
TACCGGGGTCCGAAGA-1  G12V
TATACCTCATAGGCGA-1  G12V
CTGCCATGTATTGGCT-1  G12V
AGTGTTGTCAGGAGAC-1  G12V
AATGCCAAGCAGTACG-1  G12V
TTTGTTGAGGCTTCCG-1  G12V
ACAGAAATCCATTCGC-1  G12V
GTCATTTCAACTAGAA-1  G12V
GAGTGTTGTGAGATTA-1  G12V
TCTTGCGCATAACGGG-1  G12V
ACATCGATCGGCATCG-1  G12V
CATCCCACACACAGCC-1  G12V
GACAGCCGTAGCTGTT-1  G12V
CACAGATGTTAGGCCC-1  G12V
AACGGGATCACAACCA-1  G12V

If you want to add metadata to a Seurat object the new metadata needs to have the same length as the old metadata (== number of cells).
For an example, if you have a 8 cells c("a", "b", "c", "d", "e", "f", "g") and you want to be able to highlight the 4 cells with the "barcodes" c("a", "c", "f", "g") in a umap plot you have to add a new grouping column:
# The Seurat object is called "ExampleData" and columns can be directly adressed using "$"
# The barcodes of all cells are stored as the rownames of the meta.data slot
ExampleData$cells_of_interest <- ifelse(rownames(ExampleData@meta.data) %in% c("a", "c", "f", "g"), "interesting", "boring")

From this you will get a new column in the meta data with the values:
interesing, boring, interesting, boring, boring, interesting, interesting

This way you can use this column in DimPlot:
DimPlot(ExampleData, reduction = "umap", group.by = "cells_of_interest")

With two colors in the plot one called "interesting" and one "boring"

so I would do ptx_human_patient.1$cells_of_interest <- ifelse(ptx_human_patient.1$barcode %in% c("", "", "", ""), "", "") I don't know what you mean by 8 cells vs. wanting to show 4 cells. I want to show all the barcodes from my excel file in the umap. I've uploaded the excel file by doing metadata = read.table("metadata.csv", sep=';', header=T, stringsAsFactors=F)

As I understand you want to highlight a subset of your cells using the barcodes you uploaded from file. Unsing barcode names in the original metadat you can tell Seurat wich barcodes are of interest for you and then highlight them.

yes that is what I want to do, I just don't understand how to go about it.

@mmpp I added some clarification to the answer, especially regarding where in the meta data the barcodes are stored.

I'm not sure my question is being understood correctly. I already have a file with the barcodes that I want to show on the umap. I don't want to show only four cells, I want to show all of the barcodes from my file.

I made an example with highlighting 4 out of 8 cells to make it easier to understand. You can easily apply this approach to any number of barcodes

i understand that but how do i connect my barcode file to the metadata? because right now, i'm not sure what it is considering boring or interesting.

if you replace the vector c("a", "c", "f", "g") from the example with the barcode column from your excel these barcodes will be called "interesting" and the rest "boring". Just play around with this a little and it will become clearer.

okay... i have a lot of barcodes to map, i have to type out every single one of them out? is there not an easier to do it by reading the excel file? i appreciate all your help.

You dont have to type the barcodes. Just read the file bc <- read.xlsx("Excel.xlsx") and then use the column with the barcodes (bc$barcode) in the place where I used the vector. All this is very basic R so if you are having problems you should look for an introductory tutorial

i did ptx_human_patient.1$cells_of_interest <- ifelse(rownames(ptx_human_patient.1@meta.data) %in% c(SD_52_JZ33_KRAS$barcode), "interesting", "boring") and the whole column just says "boring".

It may be that Seurat is adding a number after the barcode. In this case you need to add this number to the barcodes from the Excel file.

From my understanding, you are just trying to add some new information to the metadata of your Seurat object. To do so you can just add the column to meta.data slots in your Seurat object.
(https://github.com/satijalab/seurat/wiki/Seurat#slots)
Something like this:
myBarcode = rownames(seurat_obj@meta.data) #get barcode from seurat
mutation.df = mutation.df[match(myBarcode, mutation.df$barcode), ] #match the order of seurat barcode with you data
seurat_obj$mutation = mutation.df$muTaTionsTaTus #put mutation into metadata

I tried to add the barcodes but it didn't seem to work. I did it by ptx_human_patient.1$KRAS ="metadata.1"

Can you describe your question more specifically? Is there a barcode that you want to plot for each of the cell, or just a subset of cells? Can you post several of your scRNA barcodes here?

I copied some of the barcodes in the question. They are just a subset of cells. I have KRAS variants and would like to map on the umap. I also tried this ptx_human_patient.1 <- AddMetaData(object = ptx_human_patient.1, metadata = metadata.1, col.name = "variant")

I think @PPK's answer is correct. Before you add the new metadata column to Seurat, you need to make sure that they have the same number of rows and the cell barcodes are in the same order.

I'm not sure how to do that and I don't understand how it's done based on the response.

Basically to combine two data frames by common rows. Maybe you can try this:
seurat_obj@meta.data <- left_join(seurat_obj@meta.data, excelDf, by = "barcode"). Here excelDf is your dataframe from the excel. Then print the seurat_obj@meta.data to check if it is the desired results. I haven't used Seurat for a while. Let me know if you get any error message.

I get this error: Error: by can't contain join column barcode which is missing from LHS

this is how i upload the excel file : metadata = read.table("SD-52-JZ33_KRAS.csv", sep=';', header=T, stringsAsFactors=F)

sorry about this what is the column for barcode name in seurat_obj@meta.data? looks like it is not barcode. Can you print the data frame to confirm?

barcode            mutation    <chr>              <chr>     1 GAGATGGCACTAGGCC-1 G12V      2 CTGAGCGAGACTCAAA-1 G12V      3 TACCGGGGTCCGAAGA-1 G12V      4 TATACCTCATAGGCGA-1 G12V      5 CTGCCATGTATTGGCT-1 G12V      6 AGTGTTGTCAGGAGAC-1 G12V      7 AATGCCAAGCAGTACG-1 G12V      8 TTTGTTGAGGCTTCCG-1 G12V      9 ACAGAAATCCATTCGC-1 G12V     10 GTCATTTCAACTAGAA-1 G12V this is the dataframe of my excel file.

Can you try seurat_obj@meta.data?

dput(head(ptx_human_patient.1@meta.data)[1:3]) structure(list(orig.ident = c("human_mouse_alignment_ptx", "human_mouse_alignment_ptx",  "human_mouse_alignment_ptx", "human_mouse_alignment_ptx", "human_mouse_alignment_ptx",  "human_mouse_alignment_ptx"), nCount_RNA = c(15412, 21185, 5867,  16330, 14714, 12944), nFeature_RNA = c(3546L, 4221L, 1815L, 4807L,  4341L, 3055L)), row.names = c("G_AAACCCACAGGTCAAG", "G_AAACCCATCATCACCC",  "G_AAACGAAAGAGCAGCT", "G_AAACGAAGTCGTGGTC", "G_AAACGAATCCCGAGGT",  "G_AAACGCTAGGGAACAA"), class = "data.frame")

I did the output of dput(head(ptx_human_patient.1@meta.data)[1:3])

sorry about my previous line of code. From your output the cell barcodes are actually in the rownames of meta.data. But the barcodes look like "G_AAACGCTAGGGAACAA", which is different from your mutation data frame (like GTCATTTCAACTAGAA-1). So I think you first need to know what causes the difference and try to make them the same.

I got the mutation barcodesAs  from cellranger not seurat, that is the way they are given. As for the barcodes I have to add a letter in front of each barcode before I can merge all counts.


text:
Available phenotype data from GTEx
I am using gene expression data from the The Genotype-Tissue Expression (GTEx) project, available here: http://gtexportal.org. I would also like to use phenotype data for each sample (especially age, sex, race, body mass index, smoking history).
Access to such information is controlled by the NIH, so you need to apply for access.
Before I go through that probably lengthy process, I would like to know if it is worth it. Are any of the phenotype data age, sex, race, body mass index, smoking history available on a per-sample level, if you have access? I could not find this information anywhere.

If helpful for anyone else, I found that all and much more are available. See below 2 links.
https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/dataset.cgi?study_id=phs000424.v8.p2&phv=169091&phd=3910&pha=&pht=2742&phvf=&phdf=&phaf=&phtf=&dssp=1&consent=&temp=1
https://ftp.ncbi.nlm.nih.gov/dbgap/studies/phs000424/phs000424.v8.p2/pheno_variable_summaries/phs000424.v8.pht002742.v8.p2.GTEx_Subject_Phenotypes.var_report.xml


text:
How to run a one-time 'setup' step for amrfinder in snakemake using conda
I'm trying to add amrfinder running within conda to a snakemake workflow, working on a slurm HPC.
The first time you run amrfinder, you have to run amrfinder -u to update/download the database of genes/proteins which it searches for. Since I'm using the --use-conda flag of snakemake, the first time I run the workflow, it creates a new env for the conda package, so this command has to be run before that copy of amrfinder will work.
My question is: how do I run amrfinder -u within snakemake? If I do it as part of the standard amrfinder command, it runs as many times as there are samples, and as all the jobs are looking at the same database, when it updates for one job, it breaks the other jobs.
I can't figure out how to run this command in a 'snakemake-y' way, as it has no input, it's just literally amrfinder -u.
Of course, I can make it work with a hack/workaround, but I'm interested in a snakemake way of doing it.
This is the main snakemake rule for running amr_finder_plus, the new rule, to run amrfinder -u needs to run before this
rule amr_finder_plus:
    input:
        assembly = 'my_assembly.fasta'
    output:
        amr_finder_plus_results = 'my_amr_finder_plus_results.tsv'
    conda:
        '../../envs/amrfinderplus.yaml'
    shell:
        'amrfinder -n {input.assembly} -O Salmonella --output {output.amr_finder_plus_results} --threads 4'

There are two approaches you can take:

There's a --conda-create-envs-only option that you can use to exit after creating the envs. You should get information about where the envs are so you can run amrfinder -u.
Add a new rule called "amr_setup" or something like that that runs:

.
amrfinder -u
touch amrfinder.setup

Then put amrfinder.setup into your snakemake workflow as an input file into the amr_finder_plus rule. Either of these will work


text:
After KEGG and GO analysis, how to make tables+phylogenetic trees
hope everyone is ok. I used Trinity to do a de novo transcriptome assembly, then blastp/blastx and then used Blast2GO software to do KEGG and GO analysis. So i got some txt files with header :
for GO
Sequence name Sequence desc. Sequence length Hit desc. Hit ACC E-Value Similarity Bit-Score Alignment length Positives
for KEGG
Seq Pathway Pathway ID #Enzs in Pathway Enzyme
All well so far. Then i was asked to do 3 tasks.

Tables with terpene synthases, (monoterpenes, sesquiterpenes, diterpenes and triterpenes) and corresponding phylogenetic trees including other known terpene synthases from all categories from other plant species.
Tables with Cytochrome P450s enzymes implicated in terpene biosynthesis (categorization per clades and clans) and corresponding phylogenetic trees including other known terpene CYPs from all clades and clans from other plant species.
From the graphs corresponding to the different KEGG and GO ontology analyses, export the data for the GO analysis for all three parameters (categories, names, blasts, total number of contigs and percentage) to an excel file to be able to transform them in other formats for future publication.

I have no idea how to proceed. Any guidance? Thanks

This answer which was given for my question would be of some help regarding the phylogenetic tree given by @Michael https://bioinformatics.stackexchange.com/questions/13481/parsimony-and-maximum-likelihood-tree-comparison-in-r

I'm going to answer 1+2 because 3 is hard to understand and maybe trivial?
Each of 1+2 seems to be have these sub-tasks:

identify the contigs in your assembly that have the genes of interest.
pull out the gene sequences from the contigs (you should have intermediate tblastx output somewhere with coordinates, if you don't you can just re-blast those contigs).
possibly translate those in silico to protein sequences.
tblastx (gene) or blastp (protein) those gene/protein sequences against whatever db is relevant.
pull sequences for the hits above some threshold for each gene from the db.
perform a multiple sequence alignment for each gene against the hits.
plug the MSA into whatever phylogenetic tool you like. I used phyml a lot back in the day, don't know what the kids like now.
make a visualization of the tree, probably?


text:
Running snakemake with multiple conda environments without rebuilding environments every time
I have some snakemake pipelines that have rules which run in different environments. Is it possible to run the pipeline without having to rebuild the environments from scratch every time when run in a different directory?
This is the way I would normally run snakemake:
snakemake --directory <my_output_dir> --use-conda
My snakefile might have rules that need to use different environments:
rule A:
   ....
   conda: envA.yml
   ....

rule B:
   ....
   conda: envB.yml
   ....

So if I wanted to do two different analyses with the same pipeline, is there anyway that I can point them to the same environments without having to install everything separately in each output directory?
The snakemake documentation makes it seem like they don't currently support this. You can either use 1) --conda-prefix when running snakemake and then point to a single env or 2) use the --use-conda directive and the env.yml files by rule, in which case it will create new envs any time you run it in a new output directory. So my case doesn't seem to be covered, but does anyone have a workaround to suggest? I often find conflicts between some of the R packages I use and python packages, so I tend to keep them in their own environments to avoid headache.
I'm starting to get neck-deep in conda environments at this point which is taking up space from my disk quota, but also I am bothered by the fact that I may actually want to run analysis of data from different projects with the same environments. If the environments are rebuilt at each execution, I do not know that they are the same, short of specifying the exact version info for everything I use in the env.yaml files (and even then, they could be slightly different if the repo has changed).
I have looked at this similar question: Running Snakemake in one single conda env
in which it is suggested there is a way to do such a thing, but I'm afraid I don't understand how to implement it without an example.
Thanks

Each YAML file will always create its own environment. The packages in each aren't downloaded and installed each time, if they're already installed elsewhere they're just hard-linked in. You'll want to use both --use-conda and --conda-prefix. The latter will point to the base directory in which all conda environments will be stored. This has the benefit that running the same workflow in different places won't result in a new env being created (by default snakemake will do so inside a .snakemake directory otherwise).

Wow, thank you. I had no idea that these two flags could be used together (I thought --conda-prefix could only be used with a named pre-built single environment) -- that solves the problem. I see what you mean about the hard-linking but it still ends up looking quite messy if you have all these .snakemake directories everywhere and it also seems to slow down some conda commands (like conda env list), at least for me.

@vantom, while I am not very familiar with snakemake on a whole, I do have extensive knowledge writing workflows in WDL as well as CWL. The approach both of these languages take to overcome a similar issue that you are describing, is to use containerization through ie docker. This approach allows you to bundle the specific tools you need for a specific task, making different images for different tasks. There are several large repositories of images online ie: Docker Hub, or more biologically focused: Biocontainers
Looking through the documentation, it does appear that snakemake supports  running your rule from within the context of a docker image. Without knowing how your compute environment is setup, I am not sure if there is any additional configuration you would need, however the documentation seems to suggest that it should be as easy as describing the name of the container. As long as you can run docker, (or an alternative like singularity) you should be able to avoid the current headache you are facing.
rule NAME:
    input:
        "table.txt"
    output:
        "plots/myplot.pdf"
    containers:
        "docker://joseespinosa/docker-r-ggplot2"
    script:
        "scripts/plot-stuff.R"

Thank you for the tip! I'm not currently using Docker, but it's useful to know this for the future.

For reproducibility, you should always used a tagged version of a container, that is:
docker://joseespinosa/docker-r-ggplot2:v0.1
If you don't it will pull/build the :latest tag and you run the risk your pipeline could break on an update.
See the tags menu in DockerHub; e.g. https://hub.docker.com/r/joseespinosa/docker-r-ggplot2/tags


text:
Why is library(Seurat) causing my R session to have a fatal error?
I recently tried to load library(Seurat) after installing Seurat package. It caused my R studio to have a fatal error and abort my session. I have successfully downloaded it before. I have uninstalled R and R studio and reinstalled, restarted my computer. I am new to R too, so I'm not sure what I'm doing.
I am following the directions provided on the satijalab.org website on how to install Seurat.
Here's a screenshot of the error message.
I ran install.packages("Seurat") in R. R is also crashing as soon as I run library(Seurat). Here's a screenshot right before I hit run, and when it crashes, R application just closes.

hi @DianaatDuke, looking at the screen shot, seems like it downloaded the package..  i am surprised it doesn't try to install other dependencies, since you installed a new R or Rstudio? can you load other libraries, save library(uwot) , since this is needed for Seurat too

Yes, no error message when i run the code library(uwot).

Okay it might be Seurat is not checking for the installation of core dependencies

I did the same analysis last week and it was actually perfectly fine. So do you think something changed since then?

i would go crazy if that's the case lol... can you try opening R , preferably use a terminal , powershell whichever you are on, and try to load Seurat.. just check it's not Rstudio

I am going crazy, it's been 2.5 hours of me trying to figure out what's wrong. R is also crashing. I can run this without problem -  install.packages("Seurat"). As soon as I run library(Seurat), R just crashes.

it can be solved..ok do you always run R from 1 directory? did you try removing the .Rhistory .Rdata files?https://stat.ethz.ch/pipermail/r-sig-mac/2012-June/009326.html

I tried to remove .Rhistory in R studio on the bottom right screen of (under Files tab). but I didn't try removing .Renviron    - where can I find and delete that?

I'm using my work mac laptop. I actually tried to run this exercise on my personal Mac, and had the same problems. I can borrow a PC laptop and try to do it on there? but ultimately, I need this to run on my mac laptop.

it's ok, mac makes it easier. So open a terminal (go to applications /utilities or just search). then do ls .R* ; you should see .Rapp.history .Renviron .Rhistory or maybe something more, delete these by doing rm .Rhistory etc..

I suggest those above because from your R, i can see you are loading .Rapp.history.. whatever that might be. Once all deleted;  type R in the terminal ; and R will be loaded in a terminal. try library(Seurat) then

this is the message I get when I type in .R*
GI-C02ZC1G0:~ dl319$ .R*
-bash: .Rapp.history: command not found

Let us continue this discussion in chat.

@StupidWolf helped me figured it out! Needed to remove r-miniconda and then reinstall Seurat package.

If you haven't already done so I would make a complete installation of Anaconda, particularly good with bash. The I would either install using Anaconda Navigator, or ...
# Create an R environment and activate it
conda create -n r_env r-essentials r-base python=3.7
conda activate r_env

# install Reurat
conda install -c bioconda r-seurat

# installing standard packages
conda install -c r r-randomforest --yes

If you have a Mac the instructions are here. In some systems you need to just activate r_env .
Please note if you've installed e.g. python3 with pip a conflict occurs with the Anaconda package.

I just saw this question. The solution to your problem is here

Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.


text:
10X cellranger error during count
I successfully installed cellranger and performed a test run. Even, I installed bcl2fastq. Everything looked perfect. Thus, I ran cellranger count. I ran into an error during the cellranger count run. I put command as suggested in the tutorial as follows (I double-checked any typo that can cause error).
tutorial (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/tutorial_ct)
cellranger count --id=run_count_1kpbmcs \
--fastqs=/mnt/home/user.name/yard/run_cellranger_count/pbmc_1k_v3_fastqs \
--sample=pbmc_1k_v3 \
--transcriptome=/mnt/home/user.name/yard/run_cellranger_count/refdata-cellranger-GRCh38-3.0.0

Initially, the process went smoothly, but it suddenly found error as follows. It's not likely due to memory issues.

[error] Pipestance failed. Error log at:
run_count_1kpbmcs/SC_RNA_COUNTER_CS/SC_RNA_COUNTER/_BASIC_SC_RNA_COUNTER/ALIGN_READS/fork0/chnk00-u27879f31e3/_errors

Log message:
Traceback (most recent call last):
  File "/ufrc/zhangw/my.kim/yard/apps/cellranger-3.1.0/martian-cs/v3.2.3/adapters/python/martian_shell.py", line 589, in _main
    stage.main()
  File "/ufrc/zhangw/my.kim/yard/apps/cellranger-3.1.0/martian-cs/v3.2.3/adapters/python/martian_shell.py", line 554, in main
    self._run(lambda: self._module.main(args, outs))
  File "/ufrc/zhangw/my.kim/yard/apps/cellranger-3.1.0/martian-cs/v3.2.3/adapters/python/martian_shell.py", line 523, in _run
    cmd()
  File "/ufrc/zhangw/my.kim/yard/apps/cellranger-3.1.0/martian-cs/v3.2.3/adapters/python/martian_shell.py", line 554, in <lambda>
    self._run(lambda: self._module.main(args, outs))
  File "/ufrc/zhangw/my.kim/yard/apps/cellranger-3.1.0/cellranger-cs/3.1.0/mro/stages/counter/align_reads/__init__.py", line 128, in main
    align_reads(args, outs)
  File "/ufrc/zhangw/my.kim/yard/apps/cellranger-3.1.0/cellranger-cs/3.1.0/mro/stages/counter/align_reads/__init__.py", line 67, in align_reads
    read_group_tags=tk_bam.make_star_rg_header(args.read_group))
  File "/ufrc/zhangw/my.kim/yard/apps/cellranger-3.1.0/cellranger-cs/3.1.0/lib/python/cellranger/reference.py", line 616, in align
    raise cr_io.CRCalledProcessError(e.msg + ' Check STAR logs for errors: %s .' % star_log)
CRCalledProcessError: Process returned error code -9: STAR . Check STAR logs for errors: /ufrc/zhangw/my.kim/yard/run_cellranger_count/run_count_1kpbmcs/SC_RNA_COUNTER_CS/SC_RNA_COUNTER/_BASIC_SC_RNA_COUNTER/ALIGN_READS/fork0/chnk00-u27879f31e3/files/Log.out .

Waiting 6 seconds for UI to do final refresh.
Pipestance failed. Use --noexit option to keep UI running after failure.

2020-04-21 13:48:37 Shutting down.

Meanwhile, I checked STAR logs for errors according to Haci said, and I found these,
STAR version=                 STAR_2.5.1b
STAR compilation time,server,dir=Wed May 8 16:13:25 PDT 2019 :/mnt/home/jenkins/workspace/module-STAR-5dda596---latest/sake/modules/STAR/5dda596/source
##### DEFAULT parameters:
versionSTAR                       20201
versionGenome                     20101   20200
parametersFiles                   -
sysShell                          -
runMode                           alignReads
runThreadN                        1
runDirPerm                        User_RWX
runRNGseed                        777
genomeDir                         ./GenomeDir/
genomeLoad                        NoSharedMemory
genomeFastaFiles                  -
genomeSAindexNbases               14
genomeChrBinNbits                 18
genomeSAsparseD                   1
genomeSuffixLengthMax             18446744073709551615
readFilesIn                       Read1   Read2
readFilesCommand                  -
readMatesLengthsIn                NotEqual
readMapNumber                     18446744073709551615
readNameSeparator                 /
inputBAMfile                      -
bamRemoveDuplicatesType           -
bamRemoveDuplicatesMate2basesN    0
limitGenomeGenerateRAM            31000000000
limitIObufferSize                 150000000
limitOutSAMoneReadBytes           100000
limitOutSJcollapsed               1000000
limitOutSJoneRead                 1000
limitBAMsortRAM                   0
limitSjdbInsertNsj                1000000
outFileNamePrefix                 ./
outTmpDir                         -
outReadsUnmapped                  None
outQSconversionAdd                0
outMultimapperOrder               Old_2.4
outSAMtype                        SAM
outSAMmode                        Full
outSAMstrandField                 None
outSAMattributes                  Standard
outSAMunmapped                    None
outSAMorder                       Paired
outSAMprimaryFlag                 OneBestScore
outSAMreadID                      Standard
outSAMmapqUnique                  255
outSAMflagOR                      0
outSAMflagAND                     65535
outSAMattrRGline                  -
outSAMheaderHD                    -
outSAMheaderPG                    -
outSAMheaderCommentFile           -
outBAMcompression                 1
outBAMsortingThreadN              0
outSAMfilter                      None
outSAMmultNmax                    18446744073709551615
outSAMattrIHstart                 1
outSJfilterReads                  All
outSJfilterCountUniqueMin         3   1   1   1
outSJfilterCountTotalMin          3   1   1   1
outSJfilterOverhangMin            30   12   12   12
outSJfilterDistToOtherSJmin       10   0   5   10
outSJfilterIntronMaxVsReadN       50000   100000   200000
outWigType                        None
outWigStrand                      Stranded
outWigReferencesPrefix            -
outWigNorm                        RPM
outFilterType                     Normal
outFilterMultimapNmax             10
outFilterMultimapScoreRange       1
outFilterScoreMin                 0
outFilterScoreMinOverLread        0.66
outFilterMatchNmin                0
outFilterMatchNminOverLread       0.66
outFilterMismatchNmax             10
outFilterMismatchNoverLmax        0.3
outFilterMismatchNoverReadLmax    1
outFilterIntronMotifs             None
clip5pNbases                      0
clip3pNbases                      0
clip3pAfterAdapterNbases          0
clip3pAdapterSeq                  -
clip3pAdapterMMp                  0.1
winBinNbits                       16
winAnchorDistNbins                9
winFlankNbins                     4
winAnchorMultimapNmax             50
scoreGap                          0
scoreGapNoncan                    -8
scoreGapGCAG                      -4
scoreGapATAC                      -8
scoreStitchSJshift                1
scoreGenomicLengthLog2scale       -0.25
scoreDelBase                      -2
scoreDelOpen                      -2
scoreInsOpen                      -2
scoreInsBase                      -2
seedSearchLmax                    0
seedSearchStartLmax               50
seedSearchStartLmaxOverLread      1
seedPerReadNmax                   1000
seedPerWindowNmax                 50
seedNoneLociPerWindow             10
seedMultimapNmax                  10000
alignIntronMin                    21
alignIntronMax                    0
alignMatesGapMax                  0
alignTranscriptsPerReadNmax       10000
alignSJoverhangMin                5
alignSJDBoverhangMin              3
alignSJstitchMismatchNmax         0   -1   0   0
alignSplicedMateMapLmin           0
alignSplicedMateMapLminOverLmate    0.66
alignWindowsPerReadNmax           10000
alignTranscriptsPerWindowNmax     100
alignEndsType                     Local
alignSoftClipAtReferenceEnds      Yes
chimSegmentMin                    0
chimScoreMin                      0
chimScoreDropMax                  20
chimScoreSeparation               10
chimScoreJunctionNonGTAG          -1
chimJunctionOverhangMin           20
chimOutType                       SeparateSAMold
chimFilter                        banGenomicN
chimSegmentReadGapMax             0
sjdbFileChrStartEnd               -
sjdbGTFfile                       -
sjdbGTFchrPrefix                  -
sjdbGTFfeatureExon                exon
sjdbGTFtagExonParentTranscript    transcript_id
sjdbGTFtagExonParentGene          gene_id
sjdbOverhang                      100
sjdbScore                         2
sjdbInsertSave                    Basic
quantMode                         -
quantTranscriptomeBAMcompression    1
quantTranscriptomeBan             IndelSoftclipSingleend
twopass1readsN                    18446744073709551615
twopassMode                       None
##### Command Line:
STAR --genomeDir /ufrc/zhangw/my.kim/yard/run_cellranger_count/refdata-cellranger-GRCh38-3.0.0/star --outSAMmultNmax -1 --runThreadN 4 --readNameSeparator space --outSAMunmapped Within KeepPairs --outSAMtype SAM --outStd SAM --outSAMorder PairedKeepInputOrder --outSAMattrRGline ID:run_count_1kpbmcs:0:1:HFWFVDMXX:1 SM:run_count_1kpbmcs LB:0.1 PU:run_count_1kpbmcs:0:1:HFWFVDMXX:1 PL:ILLUMINA --readFilesIn /dev/fd/63
##### Initial USER parameters from Command Line:
outStd                            SAM
###### All USER parameters from Command Line:
genomeDir                     /ufrc/zhangw/my.kim/yard/run_cellranger_count/refdata-cellranger-GRCh38-3.0.0/star     ~RE-DEFINED
outSAMmultNmax                18446744073709551615     ~RE-DEFINED
runThreadN                    4     ~RE-DEFINED
readNameSeparator             space        ~RE-DEFINED
outSAMunmapped                Within   KeepPairs        ~RE-DEFINED
outSAMtype                    SAM        ~RE-DEFINED
outStd                        SAM     ~RE-DEFINED
outSAMorder                   PairedKeepInputOrder     ~RE-DEFINED
outSAMattrRGline              ID:run_count_1kpbmcs:0:1:HFWFVDMXX:1   SM:run_count_1kpbmcs   LB:0.1   PU:run_count_1kpbmcs:0:1:HFWFVDMXX:1   PL:ILLUMINA        ~RE-DEFINED
readFilesIn                   /dev/fd/63        ~RE-DEFINED
##### Finished reading parameters from all sources

##### Final user re-defined parameters-----------------:
runThreadN                        4
genomeDir                         /ufrc/zhangw/my.kim/yard/run_cellranger_count/refdata-cellranger-GRCh38-3.0.0/star
readFilesIn                       /dev/fd/63
readNameSeparator                 space
outStd                            SAM
outSAMtype                        SAM
outSAMunmapped                    Within   KeepPairs
outSAMorder                       PairedKeepInputOrder
outSAMattrRGline                  ID:run_count_1kpbmcs:0:1:HFWFVDMXX:1   SM:run_count_1kpbmcs   LB:0.1   PU:run_count_1kpbmcs:0:1:HFWFVDMXX:1   PL:ILLUMINA
outSAMmultNmax                    18446744073709551615

-------------------------------
##### Final effective command line:
STAR   --runThreadN 4   --genomeDir /ufrc/zhangw/my.kim/yard/run_cellranger_count/refdata-cellranger-GRCh38-3.0.0/star   --readFilesIn /dev/fd/63      --readNameSeparator space      --outStd SAM   --outSAMtype SAM      --outSAMunmapped Within   KeepPairs      --outSAMorder PairedKeepInputOrder   --outSAMattrRGline ID:run_count_1kpbmcs:0:1:HFWFVDMXX:1   SM:run_count_1kpbmcs   LB:0.1   PU:run_count_1kpbmcs:0:1:HFWFVDMXX:1   PL:ILLUMINA      --outSAMmultNmax 18446744073709551615

##### Final parameters after user input--------------------------------:
versionSTAR                       20201
versionGenome                     20101   20200
parametersFiles                   -
sysShell                          -
runMode                           alignReads
runThreadN                        4
runDirPerm                        User_RWX
runRNGseed                        777
genomeDir                         /ufrc/zhangw/my.kim/yard/run_cellranger_count/refdata-cellranger-GRCh38-3.0.0/star
genomeLoad                        NoSharedMemory
genomeFastaFiles                  -
genomeSAindexNbases               14
genomeChrBinNbits                 18
genomeSAsparseD                   1
genomeSuffixLengthMax             18446744073709551615
readFilesIn                       /dev/fd/63
readFilesCommand                  -
readMatesLengthsIn                NotEqual
readMapNumber                     18446744073709551615
readNameSeparator                 space
inputBAMfile                      -
bamRemoveDuplicatesType           -
bamRemoveDuplicatesMate2basesN    0
limitGenomeGenerateRAM            31000000000
limitIObufferSize                 150000000
limitOutSAMoneReadBytes           100000
limitOutSJcollapsed               1000000
limitOutSJoneRead                 1000
limitBAMsortRAM                   0
limitSjdbInsertNsj                1000000
outFileNamePrefix                 ./
outTmpDir                         -
outStd                            SAM
outReadsUnmapped                  None
outQSconversionAdd                0
outMultimapperOrder               Old_2.4
outSAMtype                        SAM
outSAMmode                        Full
outSAMstrandField                 None
outSAMattributes                  Standard
outSAMunmapped                    Within   KeepPairs
outSAMorder                       PairedKeepInputOrder
outSAMprimaryFlag                 OneBestScore
outSAMreadID                      Standard
outSAMmapqUnique                  255
outSAMflagOR                      0
outSAMflagAND                     65535
outSAMattrRGline                  ID:run_count_1kpbmcs:0:1:HFWFVDMXX:1   SM:run_count_1kpbmcs   LB:0.1   PU:run_count_1kpbmcs:0:1:HFWFVDMXX:1   PL:ILLUMINA
outSAMheaderHD                    -
outSAMheaderPG                    -
outSAMheaderCommentFile           -
outBAMcompression                 1
outBAMsortingThreadN              0
outSAMfilter                      None
outSAMmultNmax                    18446744073709551615
outSAMattrIHstart                 1
outSJfilterReads                  All
outSJfilterCountUniqueMin         3   1   1   1
outSJfilterCountTotalMin          3   1   1   1
outSJfilterOverhangMin            30   12   12   12
outSJfilterDistToOtherSJmin       10   0   5   10
outSJfilterIntronMaxVsReadN       50000   100000   200000
outWigType                        None
outWigStrand                      Stranded
outWigReferencesPrefix            -
outWigNorm                        RPM
outFilterType                     Normal
outFilterMultimapNmax             10
outFilterMultimapScoreRange       1
outFilterScoreMin                 0
outFilterScoreMinOverLread        0.66
outFilterMatchNmin                0
outFilterMatchNminOverLread       0.66
outFilterMismatchNmax             10
outFilterMismatchNoverLmax        0.3
outFilterMismatchNoverReadLmax    1
outFilterIntronMotifs             None
clip5pNbases                      0
clip3pNbases                      0
clip3pAfterAdapterNbases          0
clip3pAdapterSeq                  -
clip3pAdapterMMp                  0.1
winBinNbits                       16
winAnchorDistNbins                9
winFlankNbins                     4
winAnchorMultimapNmax             50
scoreGap                          0
scoreGapNoncan                    -8
scoreGapGCAG                      -4
scoreGapATAC                      -8
scoreStitchSJshift                1
scoreGenomicLengthLog2scale       -0.25
scoreDelBase                      -2
scoreDelOpen                      -2
scoreInsOpen                      -2
scoreInsBase                      -2
seedSearchLmax                    0
seedSearchStartLmax               50
seedSearchStartLmaxOverLread      1
seedPerReadNmax                   1000
seedPerWindowNmax                 50
seedNoneLociPerWindow             10
seedMultimapNmax                  10000
alignIntronMin                    21
alignIntronMax                    0
alignMatesGapMax                  0
alignTranscriptsPerReadNmax       10000
alignSJoverhangMin                5
alignSJDBoverhangMin              3
alignSJstitchMismatchNmax         0   -1   0   0
alignSplicedMateMapLmin           0
alignSplicedMateMapLminOverLmate    0.66
alignWindowsPerReadNmax           10000
alignTranscriptsPerWindowNmax     100
alignEndsType                     Local
alignSoftClipAtReferenceEnds      Yes
chimSegmentMin                    0
chimScoreMin                      0
chimScoreDropMax                  20
chimScoreSeparation               10
chimScoreJunctionNonGTAG          -1
chimJunctionOverhangMin           20
chimOutType                       SeparateSAMold
chimFilter                        banGenomicN
chimSegmentReadGapMax             0
sjdbFileChrStartEnd               -
sjdbGTFfile                       -
sjdbGTFchrPrefix                  -
sjdbGTFfeatureExon                exon
sjdbGTFtagExonParentTranscript    transcript_id
sjdbGTFtagExonParentGene          gene_id
sjdbOverhang                      100
sjdbScore                         2
sjdbInsertSave                    Basic
quantMode                         -
quantTranscriptomeBAMcompression    1
quantTranscriptomeBan             IndelSoftclipSingleend
twopass1readsN                    18446744073709551615
twopassMode                       None
----------------------------------------

WARNING --outSAMattrRG defines a read group, therefore STAR will output RG attribute
Finished loading and checking parameters
Reading genome generation parameters:
versionGenome                 20201        ~RE-DEFINED
genomeFastaFiles              /mnt/yard2/pat/cr/new-refs/reference_buildscripts2/GRCh38/fasta/genome.fa        ~RE-DEFINED
genomeSAindexNbases           14     ~RE-DEFINED
genomeChrBinNbits             18     ~RE-DEFINED
genomeSAsparseD               3     ~RE-DEFINED
sjdbOverhang                  100     ~RE-DEFINED
sjdbFileChrStartEnd           -        ~RE-DEFINED
sjdbGTFfile                   /mnt/yard2/pat/cr/new-refs/reference_buildscripts2/GRCh38/genes/genes.gtf     ~RE-DEFINED
sjdbGTFchrPrefix              -     ~RE-DEFINED
sjdbGTFfeatureExon            exon     ~RE-DEFINED
sjdbGTFtagExonParentTranscripttranscript_id     ~RE-DEFINED
sjdbGTFtagExonParentGene      gene_id     ~RE-DEFINED
sjdbInsertSave                Basic     ~RE-DEFINED
Genome version is compatible with current STAR version
Number of real (reference) chromosmes= 194
1   1   248956422   0
2   10  133797422   249036800
3   11  135086622   382992384
4   12  133275309   518258688
5   13  114364328   651689984
6   14  107043718   766246912
7   15  101991189   873463808
8   16  90338345    975699968
9   17  83257441    1066139648
10  18  80373285    1149501440
11  19  58617616    1229979648
12  2   242193529   1288699904
13  20  64444167    1530920960
14  21  46709983    1595408384
15  22  50818468    1642332160
16  3   198295559   1693188096
17  4   190214555   1891631104
18  5   181538259   2081947648
19  6   170805979   2263613440
20  7   159345973   2434531328
21  8   145138636   2593914880
22  9   138394717   2739142656
23  MT  16569   2877554688
24  X   156040895   2877816832
25  Y   57227415    3034054656
26  KI270728.1  1872759 3091464192
27  KI270727.1  448248  3093561344
28  KI270442.1  392061  3094085632
29  KI270729.1  280839  3094609920
30  GL000225.1  211173  3095134208
31  KI270743.1  210658  3095396352
32  GL000008.2  209709  3095658496
33  GL000009.2  201709  3095920640
34  KI270747.1  198735  3096182784
35  KI270722.1  194050  3096444928
36  GL000194.1  191469  3096707072
37  KI270742.1  186739  3096969216
38  GL000205.2  185591  3097231360
39  GL000195.1  182896  3097493504
40  KI270736.1  181920  3097755648
41  KI270733.1  179772  3098017792
42  GL000224.1  179693  3098279936
43  GL000219.1  179198  3098542080
44  KI270719.1  176845  3098804224
45  GL000216.2  176608  3099066368
46  KI270712.1  176043  3099328512
47  KI270706.1  175055  3099590656
48  KI270725.1  172810  3099852800
49  KI270744.1  168472  3100114944
50  KI270734.1  165050  3100377088
51  GL000213.1  164239  3100639232
52  GL000220.1  161802  3100901376
53  KI270715.1  161471  3101163520
54  GL000218.1  161147  3101425664
55  KI270749.1  158759  3101687808
56  KI270741.1  157432  3101949952
57  GL000221.1  155397  3102212096
58  KI270716.1  153799  3102474240
59  KI270731.1  150754  3102736384
60  KI270751.1  150742  3102998528
61  KI270750.1  148850  3103260672
62  KI270519.1  138126  3103522816
63  GL000214.1  137718  3103784960
64  KI270708.1  127682  3104047104
65  KI270730.1  112551  3104309248
66  KI270438.1  112505  3104571392
67  KI270737.1  103838  3104833536
68  KI270721.1  100316  3105095680
69  KI270738.1  99375   3105357824
70  KI270748.1  93321   3105619968
71  KI270435.1  92983   3105882112
72  GL000208.1  92689   3106144256
73  KI270538.1  91309   3106406400
74  KI270756.1  79590   3106668544
75  KI270739.1  73985   3106930688
76  KI270757.1  71251   3107192832
77  KI270709.1  66860   3107454976
78  KI270746.1  66486   3107717120
79  KI270753.1  62944   3107979264
80  KI270589.1  44474   3108241408
81  KI270726.1  43739   3108503552
82  KI270735.1  42811   3108765696
83  KI270711.1  42210   3109027840
84  KI270745.1  41891   3109289984
85  KI270714.1  41717   3109552128
86  KI270732.1  41543   3109814272
87  KI270713.1  40745   3110076416
88  KI270754.1  40191   3110338560
89  KI270710.1  40176   3110600704
90  KI270717.1  40062   3110862848
91  KI270724.1  39555   3111124992
92  KI270720.1  39050   3111387136
93  KI270723.1  38115   3111649280
94  KI270718.1  38054   3111911424
95  KI270317.1  37690   3112173568
96  KI270740.1  37240   3112435712
97  KI270755.1  36723   3112697856
98  KI270707.1  32032   3112960000
99  KI270579.1  31033   3113222144
100 KI270752.1  27745   3113484288
101 KI270512.1  22689   3113746432
102 KI270322.1  21476   3114008576
103 GL000226.1  15008   3114270720
104 KI270311.1  12399   3114532864
105 KI270366.1  8320    3114795008
106 KI270511.1  8127    3115057152
107 KI270448.1  7992    3115319296
108 KI270521.1  7642    3115581440
109 KI270581.1  7046    3115843584
110 KI270582.1  6504    3116105728
111 KI270515.1  6361    3116367872
112 KI270588.1  6158    3116630016
113 KI270591.1  5796    3116892160
114 KI270522.1  5674    3117154304
115 KI270507.1  5353    3117416448
116 KI270590.1  4685    3117678592
117 KI270584.1  4513    3117940736
118 KI270320.1  4416    3118202880
119 KI270382.1  4215    3118465024
120 KI270468.1  4055    3118727168
121 KI270467.1  3920    3118989312
122 KI270362.1  3530    3119251456
123 KI270517.1  3253    3119513600
124 KI270593.1  3041    3119775744
125 KI270528.1  2983    3120037888
126 KI270587.1  2969    3120300032
127 KI270364.1  2855    3120562176
128 KI270371.1  2805    3120824320
129 KI270333.1  2699    3121086464
130 KI270374.1  2656    3121348608
131 KI270411.1  2646    3121610752
132 KI270414.1  2489    3121872896
133 KI270510.1  2415    3122135040
134 KI270390.1  2387    3122397184
135 KI270375.1  2378    3122659328
136 KI270420.1  2321    3122921472
137 KI270509.1  2318    3123183616
138 KI270315.1  2276    3123445760
139 KI270302.1  2274    3123707904
140 KI270518.1  2186    3123970048
141 KI270530.1  2168    3124232192
142 KI270304.1  2165    3124494336
143 KI270418.1  2145    3124756480
144 KI270424.1  2140    3125018624
145 KI270417.1  2043    3125280768
146 KI270508.1  1951    3125542912
147 KI270303.1  1942    3125805056
148 KI270381.1  1930    3126067200
149 KI270529.1  1899    3126329344
150 KI270425.1  1884    3126591488
151 KI270396.1  1880    3126853632
152 KI270363.1  1803    3127115776
153 KI270386.1  1788    3127377920
154 KI270465.1  1774    3127640064
155 KI270383.1  1750    3127902208
156 KI270384.1  1658    3128164352
157 KI270330.1  1652    3128426496
158 KI270372.1  1650    3128688640
159 KI270548.1  1599    3128950784
160 KI270580.1  1553    3129212928
161 KI270387.1  1537    3129475072
162 KI270391.1  1484    3129737216
163 KI270305.1  1472    3129999360
164 KI270373.1  1451    3130261504
165 KI270422.1  1445    3130523648
166 KI270316.1  1444    3130785792
167 KI270340.1  1428    3131047936
168 KI270338.1  1428    3131310080
169 KI270583.1  1400    3131572224
170 KI270334.1  1368    3131834368
171 KI270429.1  1361    3132096512
172 KI270393.1  1308    3132358656
173 KI270516.1  1300    3132620800
174 KI270389.1  1298    3132882944
175 KI270466.1  1233    3133145088
176 KI270388.1  1216    3133407232
177 KI270544.1  1202    3133669376
178 KI270310.1  1201    3133931520
179 KI270412.1  1179    3134193664
180 KI270395.1  1143    3134455808
181 KI270376.1  1136    3134717952
182 KI270337.1  1121    3134980096
183 KI270335.1  1048    3135242240
184 KI270378.1  1048    3135504384
185 KI270379.1  1045    3135766528
186 KI270329.1  1040    3136028672
187 KI270419.1  1029    3136290816
188 KI270336.1  1026    3136552960
189 KI270312.1  998 3136815104
190 KI270539.1  993 3137077248
191 KI270385.1  990 3137339392
192 KI270423.1  981 3137601536
193 KI270392.1  971 3137863680
194 KI270394.1  970 3138125824
--sjdbOverhang = 100 taken from the generated genome
Started loading the genome: Tue Apr 21 13:48:20 2020

checking Genome sizefile size: 3205385489 bytes; state: good=1 eof=0 fail=0 bad=0
checking SA sizefile size: 8651064546 bytes; state: good=1 eof=0 fail=0 bad=0
checking /SAindex sizefile size: 1565873619 bytes; state: good=1 eof=0 fail=0 bad=0
Read from SAindex: genomeSAindexNbases=14  nSAi=357913940
nGenome=3205385489;  nSAbyte=8651064546
GstrandBit=32   SA number of indices=2097227768
Shared memory is not used for genomes. Allocated a private copy of the genome.
Genome file size: 3205385489 bytes; state: good=1 eof=0 fail=0 bad=0
Loading Genome ... cat: .: 디렉터리입니다
-----------------

I googled it but it failed to solve the issue. What should I do next? Do you have any idea?
Thank you for helping me, toddler.
Additional information:

the issue seemed to be related to annotation, but MD5 for both fastq and reference file is correct.

Here are the errors related to my problem:
 (1) WARNING --outSAMattrRG defines a read group, therefore STAR will output RG attribute
 Finished loading and checking parameters
 Shared memory is not used for genomes

I'm having the same problem.
How do you solve it?
Thanks in advance

Yes. It's a memory issue so that I increased the memory as suggested by 10X genomics.

If you've answered your own question, please add this as an accepted answer, especially given that so a few other people have the same issue. This makes it easier for search engines to find the answer, and stops the bot from recycling the question to fish for more answers.

Sorry to bother you and thanks for the comment. How to complete this answer section then?


text:
What column and row naming requirements exist with Seurat (context: when loading SPLiT-Seq data)
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
[1] "AAACATCGAAACATCG_0" "AAACATCGAAACATCG_1"

$AAACATCGAACGTGAT
[1] "AAACATCGAACGTGAT_0"

$AACGTGATAAACATCG
[1] "AACGTGATAAACATCG_0" "AACGTGATAAACATCG_1"

$AACGTGATAACGTGAT

[1] "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.
UPDATE
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
[1] "AAACATCGAAACATCG.0" "AAACATCGAAACATCG.1" "AAACATCGAACGTGAT.0" "AACGTGATAAACATCG.0" "AACGTGATAAACATCG.1" "AACGTGATAACGTGAT.0"
[7] "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.

hepcat72 you have thanks @haci for their answer. Please do consider "accepting" the answer, because it appears close to what you are requesting and/or upvoting.

Read10X() works on mtx files, that is not the root of your problem. The problem is that you are pushing non-10x-data to look and work like 10x-data.
You can just use CreateSeuratObject() and feed your count matrix to this function. I am not sure if this particular function accepts mtx files, could not find this in the documentation, but in the worst case you will just need to supply a plain count matrix, count matrix being genes as rows and cells as columns.
And that could be a good idea to show the first few rows and columns of your SPLIT-Seq output so that we can have and idea.
Regarding gene identifiers, I would stick with gene symbols, to help you in cluster annotation but also using the "MT" prefix for mitochondrial genes for example.
It looks like you are providing somewhat obscure identifiers as gene identifiers, a mix of genome build, ensemble gene id and gene symbol and apparently Seurat does not appreciate the ,. While I cannot give you a definitive recipe on what can and cannot be set as gene identifiers when using Seurat, I can tell that I have worked both with gene symbols and Ensembl identifiers, separately, and both was fine with Seurat.

Yes, I converted my data to tsv where rows are genes and columns are cells. Setting aside my concerns about cell column names, don’t the row names have to be unique, especially when the genes are a collection from 2 species? I also noted that the hg38 Ensembl gene symbols already have MT- in them, but mm10 doesn’t, so I would have to address that.

In theory the human and mouse gene symbols would not be the same for the "same" (actually orthologous genes); human gene symbols would be all uppercase and those of mouse would start with an uppercase letter followed by lowercase. Anyway, you can alter gene names anyway you like, doing that before reading the data into Seurat would be the "safest" I believe.

And yes, I’m feeding CreateSeuratObject the dge_data I described above.

I also did end up modifying the column names (the cell barcode IDs) to replace the underscore with a dash. When I do that and call CellsByIdentities, it returns a single list of 7 cell barcodes instead of 3 pairs and 1 single barcode, so it seems like not having underscores in those column names is critical? I tried to address all of the row/column name issues I encountered and after doing so, the plots all look much better. For example, the violin plots now all have curvy brown columns that weren’t there before I changed the row/col names for commas and underscores.

I have edited my answer, basically I would not have commas in my gene identifiers. Secondly, I have underscores in my cell names and Seurat is fine with it, in fact I believe it was one of the Seurat functions appending underscores to cell names while merging data.

Thanks so much for your help. I will add a sample of the data to the question in a bit.

OK, I updated the question to add the requested sample.  Note, the core facility requested the cell barcodes to be used as the column names. Also, WRT hg38's and mm10's gene symbols' uniqueness... My background is more prokaryotic. Gene symbols in that realm are definitely not unique and my tentative understanding was that symbols for eukaryotes were not unique either, particularly when it comes to paralogues. However, I suspect that only unique mappings make it through the process, which mitigates uniqueness issues, but may not eliminate it totally. Does Seurat group by common gene name?

With the recent additions to your question, now it looks like it would be the best to address to the Seurat developers. They are quite active on GitHub.

Thanks. Yes, I posted a couple days ago (https://github.com/satijalab/seurat/issues/2887).  I guess I am just anxious to hear back. Your feedback has been valuable though because I wasn't sure whether I was doing something wrong or not. And it sounds like I am doing things roughly correctly and that these issues may not necessarily be due to me doing something wrong (though I'm still not entirely sure).


text:
Is there a tool for synonymous recoding of DNA sequences?
I've got a DNA sequence that I'd like to make synonymous mutations throughout, thereby preserving the amino acid sequence. Does anyone know of a tool to achieve this

I don't but it wouldn't it be easier to write one?  The other issue is that context would be useful, it is fairly common to consider synonymous mutations in fitness, particularly cancers and it would make sense to construct a null distribution.

There seems to be a bioconductor package named synmut which does just that.
Apparently, it can take codon usage into account when generating synonymous mutations.
Copy pasting an example that matches your use case from the official doc:
library("SynMut")
filepath.fasta <- system.file("extdata", "example.fasta", package = "SynMut")
filepath.csv <- system.file("extdata", "target_regions.csv", package = "SynMut")
region <- read.csv(filepath.csv)
rgd.seq <- input_seq(filepath.fasta, region)
mut.seq <- codon_random(rgd.seq, keep = TRUE)
Biostrings::writeXStringSet(get_dna(mut.seq), "rgd.fasta")

Brilliant, thanks, that package works well

DNA Chisel can introduce synonymous mutations into your sequence.
Annotate the coding sequence with the @cds constraint and the ~change objective.
Additionally, you may also optimize the codon usage to match your species with ~match_codon_usage or one of the other options.
A web interface is also provided at EGF CUBA (Collection of Useful Biological Apps): Sculpt a Sequence
Disclaimer: I'm the current maintainer of DNA Chisel


text:
What is the perferred method of optimization or energy minimization of small molecules downloaded from PubChem?
I will be docking a small library of molecules that I have downloaded from PubChem using AutoDock. I thought why not minimize their energy before docking using ChemDraw. I am thinking; what is the best method of optimization/ energy minimization of small molecules. MM2, MMF94 or other? And why?

The difference between different force-fields is not going to be major, it is the side steps which are.
When
If you are starting from a SMILES string, optimisation is a must obviously.
If you are using a 3D conformer from PubChem or even an actual sub-1 Å crystal structure from CSD, optimisation is nice for consistency.
Which
MMF94 is a solid choice. RDkit offers MMF94 (Merck force field 1994) and UFF (universal force field). The latter has some extra features, such as not kicking the bucket at atoms that aren't your usual suspects. ChemDraw offers MM2 and MMF94 IIRC, the former is very old. g16 in Gaussian (paid software) uses Gaussian Electrostatic Model, which is polarisable and will solve partial charges in a nice way (which I'll get to in a minute). OpenBabel uses MMF94 under the hood. You can even optimise with macromolecular forcefields, such as AMBER and CHARMM.
Given that the AutoDock uses a flexible ligand it should not matter too much anyway.
What else?
Autodock does not require conformers by virtue of the flexible ligand business (the sampler does rototorsions). Other programs might if they use rigid ligands (e.g. Fred) or they use a mixed approach and use a conformer library for speed (e.g. Rosetta ligand_dock). Conformer diversity is actually a bigger issue as conformers that are only slightly contorted, often even less than kB (0.6 kcal/mol), may not be sampled.
The thing that messes things is incorrect protonation. OpenBabel can protonate at a given pH, RDKit cannot. Not an issue for ligands with drug screens (cf. Lipinski rule), but is for natural substrates. However, for the protein it is essential to make sure acidic and basic residue are correctly protonated and that histidine have the correct tautomer —there are tools, c.f. MOE, that can do it but by far the best approach is reading up on the mechanism.
The next thing is dodgy Gasteiger Marsili partial charges —for Autodock 4 (and most other full/hybrid physics-based docking programs), partial charges are important: after all PDBQT has a Q in it for charge! For prepare_ligand4.py you give it a mol2 file. A major consideration is what format you are using for your molecules. This sounds a very silly thing to point out, but you'd be surprised how many trip up. An mol file (or an sdf file, which is the same but multi-entry) does not have an official partial charge column (although they are often added as comments after M END). So whatever force-field you chose make sure that the program you use gave you partial charges are not formal charges. This does not quite apply to AutoDock Vina, which is charge independent and more reliant on machine learning derived factors.
Verdict
Open babel is nice in that it can start with a SMILES, handles multiple formats including mol2 and PDBQT, can protonate at pH 7, uses MMF94 and does Gasteiger Marsili charges (and its confab algorithm is very good, not need here though).


text:
Phasing problem using "samtools phase"
so I'm having a problem using samtools phase. I'm trying to phase my .bam archives that where generated by aligning my samples with my reference genome. Following the documentation I tried the command:
samtools phase -A  my_referenceGenome.fasta   my_sample_toBePhased.bam
The -A flag stands for dropping reads with ambiguous phase. The error I get is: [W::sam_read1] Parse error at line 1. I generated .bam files with and without the SAM header but the error persists.
I feel like I am missing something basic but that isn't clearly stated anywere. Does the my_referenceGenome.fasta archive needs to be previously phased? If so, how can I do it?
Thank you so much!

This is what I see in my samtools phase (v1.10)
Usage:   samtools phase [options] <in.bam>

samtools phase only takes a bam file as an input. But you are specifying a .fasta file as an input. If you want to specify the reference genome, you should be using the argument "--reference" as described here
      --reference FILE
           Reference sequence FASTA FILE [null]

Oh, you're right! I also tried that but it didn't work. I managed to make it work and I'll add an answer shortly. But before that I want to ask a little question:  in my pipeline I also remove pcr duplicates with the fixmate/markdup steps. Is it better to phase the reads before or after removing the duplicates?

I would recommend removing duplicates first.


text:
how to use list of gene id to get cds sequence(cds fasta file have many annotation, only gene id: is same to query id)
i have a question when i want to extract cds sequence using gene id. but cds file is not just start with >gene is, it has many other annotation. the only same  is star with gene:
cds fasta:
>Zm002 cds gene:Zm1d035916 gene_biotype:protein_coding
ATCGGCAT
>Zm001 cds RefGen_v4:9:153880862:153883850:-1 gene:Zm1d048 gene_biotype:protein_coding
ATGCGGCA

gene_list
Zm1d035916
Zm1d048

how to get result like
>Zm1d035916
ATCGGCAT
>Zm1d048
ATGCGGCA

If you can use awk, you could use this simple script:
awk '
BEGIN {p=0}
NR==FNR {gid[$1];next}
/^>/{
  p=0;
  for(g in gid){
    re="\\sgene:"g"\\s"
    if($0~re){p=1;break}
  }
}p' gene_list cds.fasta

You can also put the code inside a script file:
extract_genes.awk
#!/usr/bin/awk
BEGIN {
  p = 0; # p indicates whether to print a line or not
}
# while scanning the first file, store gene ids in an array
NR==FNR {
  gid[$1];
  next
}
# check FASTA headers
/^>/ {
  p = 0;
  for(g in gid) {
    # construct regular expression to match gene id in FASTA header
    re = "\\sgene:"g"\\s";
    if ($0~re) { # a match was found!
      p = 1;
      break;
    }
  }
}
p # print line or not?

And call it like this:
awk -f extract_genes.awk gene_list cds.fasta

hi, thank you for your answer. Where to find my output file? After I run it, I don't know where to find the output file when finished. Thank you!

Happy to help. The output is written to stdout. You can forward it into a file like this: awk -f extract_genes.awk gene_list cds.fasta > output.fasta

Thank you. However, I don't know why the output.fasta file is empty. I use command like this: awk -f gene_cds_extract.awk gene_list_red_green_cpsf100.txt all_species_longset_cds.fa > example_cds.fa

That's odd. Are the input files formatted the same as in your question? And are there overlaps between the gene names and the fasta headers?

Thank you for your quick reply. input cds.fast format like this: >PNS24249 cds supercontig:Brachypodium_distachyon_v3.0:Bd1_centromere_containing_Bradi1g41430:24784:29531:-1 gene:BRADI_41430s00200v3 gene_biotype:protein_coding transcript_biotype:protein_coding description:hypothetical protein
ATGCGGCCGCGCGGCAGAGGAGACGACGGGGACGACGACGAAGACGAAGAGCCCTCGCGT
CAACGGAGCGCGCCGGAAGCCGACCCTGATCCCGAGCCCGCGCCAGCGCAGCAGCCCACA
CGAGCGCCCCTCAGCAGCCTAGTCGTCAGGCCGCCGCCGCCGCAGGACAACGGCGCCTCC     input list.txt: BRADI_41430s00200v3 like this. no gene overlap, I extract longest cds in cds.fa file

I just made a test with your input and the script above and it outputs the sequence. Which operating system do you use, and which awk version (I use Lubuntu 18.04 and GNU Awk 4.1.4)?

One thing you could try is to use a simpler (less specific) regular expression. If you replace re = "\\sgene:"g"\\s" with re = "gene:"g does it output anything?

I use mac terminal. I don't know it already have awk or not. Sorry, It's my first time to use awk in terminal.

Thank you, I try it with simpler regular expression, and it got results. Great!!! Thank you very much.

Aha! Great that it worked. Maybe awk on Mac understands special characters in regular expressions (in this case "\s" for a space) differently than awk on Linux.

Anyway, it saves a lot time to extract sequence, not just find, copy and paste, you make it more interesting. Thank you!

If you read it in python, they would take care of that:
from Bio import SeqIO
records = list(SeqIO.parse("test.fa","fasta"))

records

[SeqRecord(seq=Seq('ATCGGCAT', SingleLetterAlphabet()), id='Zm002', name='Zm002', description='Zm002 cds gene:Zm1d035916 gene_biotype:protein_coding', dbxrefs=[]),
 SeqRecord(seq=Seq('ATGCGGCA', SingleLetterAlphabet()), id='Zm001', name='Zm001', description='Zm001 cds RefGen_v4:9:153880862:153883850:-1 gene:Zm1d048 gene_biotype:protein_coding', dbxrefs=[])]

Get the gene and sequence.. maybe there's a easier to do this in biopython, this is like my 10th time using it:
def format(rec):
    des = rec.description.split()
    gene = [entry[5:] for entry in des if "gene:" in entry][0]
    return(">"+gene+"\n"+str(rec.seq))

[format(rec) for rec in records]

['>Zm1d035916\nATCGGCAT', '>Zm1d048\nATGCGGCA']

Write it out:
f=open('example.fasta','w')
f.writelines("\n".join([format(rec) for rec in records]))
f.close()

list(SeqIO.parse("example.fasta","fasta"))
[SeqRecord(seq=Seq('ATCGGCAT', SingleLetterAlphabet()), id='Zm1d035916', name='Zm1d035916', description='Zm1d035916', dbxrefs=[]),
 SeqRecord(seq=Seq('ATGCGGCA', SingleLetterAlphabet()), id='Zm1d048', name='Zm1d048', description='Zm1d048', dbxrefs=[])]

This is a very neat comprehension

lol OP chose the mighty awk =p . Thank you for the kind comment

The OP didn't appreciate the code. This comprehension is cool and pythonic f.writelines("\n".join([format(rec) for rec in records]))

Thanks! I try it and it is also work.


text:
How to downsample some of the samples in RNA-seq data?
I have 40 samples and these are into two groups. I would like to perform a differential analysis between two groups. The library size of the samples is very low. But there are two samples (GroupA_12 and GroupB_2) in which the number of counts is too large compared to other samples. Please check the library sizes below:

I wanted to downsample these two cases and then would like to perform differential analysis. May I know how to downsample those two samples.
Any help is appreciated. thanq.

You have a few options:

Downsample the fastq files and rerun the entire analysis. You can do this with seqtk sample.
Downsample the BAM files, which you can do with samtools view -s.
Divide all of the counts in the counts files by some factor and round that to an integer.

I personally prefer option 2, since it's quick and doesn't usually have any side-effects unless you use two-pass alignment or did the quantification with something like salmon. Option 3 is always tempting, but I have to say I've never been convinced that this consistently produces similar enough results to options 1 or 2.

thanks for the reply @Devon Ryan, Yes, I feel that the 3rd option is much easier. To use the 2nd option, I have to mention that the raw read counts of the transcripts were extracted from the stringtie step using prepDE.py python script [https://ccb.jhu.edu/software/stringtie/index.shtml?t=manual] Will there be any problem?

I don't use stringTie, so I can't comment with certainty.

ok. And if I have to use the 2nd option, Is this right? samtools view -s number GroupA_12.bam >  GroupA_12_new.bam What I should give the fixed number for downsampling?

If you need 20% then 1234.2. 1234 is just a seed for reproducibility, .2 is the fraction to keep.

Small question. How to decide the fraction of reads to keep? I mean 0.2 or 0.5 or something else? On what basis I should give the fraction?

You’re trying to make the library sizes similar.

yes...so any cutoff is fine?

Just eye ball a fraction from that graph.


text:
Assemble reads for a specific gene
I have a lot of unassembled sequencing data and I want to build phylogeny for some genes from these data. I can retrieved the sequene by mapping reads to a reference sequence but if the two species are too divergent, reads will not match and I won't get the full CDS. Is there a better way to deal with that? Maybe I can do that to retrieve the well conserved part and then extend this region with reads that overlap (as if it was gap-filling)?
What do you think of this approach? Do you know a software for such extension process?
Thanks

That's a good question. I remember seeing a poster about making trees directly from reads, but unfortunately, I can not find any other record than this tweet. Perhaps you can try to email David, if he would be keen to share the implementation.
While searching for any information about David's work, I run into a GitHub repository called reads2trees. It's quite old and it's not super popular, so hard to say if it can be useful. But the name matches the problem description quite well.
I also remember there were some attempts (An example paper) to get phylogeny directly from kmers, but I have never tried that on my own. I also think these approached were focused on microbial phylogienies. So again, not sure if it's going to be useful for you.
That's all I got.

Thank you very much for your suggestions! I'm going to check that, however the ideal for me would be to retrieve the sequences not just the phylogenetic tree, for example to do selection analysis.

The "reads2trees" techniques should be able to give you sequences too (it's just a fancy name for a special case of reference-based assembly anyway). The kmer-based methods certainly won't. Good luck


text:
How to map gene names to ensembl Ids in R
I am trying to convert a list of gene names into Ensembl IDs but I'm having a bit of trouble.
I've tried using the EnsDb.Hsapiens.v86 package but I always get LRG-ID values returned as well, even when using the GeneIdFilter "startsWith" thing.
My code looks like this:
targets <- as.character(targetCSV[, 2])

targetEnsemble <- ensembldb::select(
      EnsDb.Hsapiens.v86,
      filter = GeneIdFilter("ENS", "startsWith"),
      keys = targets,
      keytype = "SYMBOL",
      columns = c("SYMBOL","GENEID"))

And for example I get back:
    SYMBOL          GENEID
1     CYR61 ENSG00000142871
2      CTGF ENSG00000118523
3    AMOTL2 ENSG00000114019
4    ANKRD1 ENSG00000148677
5    ANKRD1         LRG_379
6    IGFBP3 ENSG00000146674
7        F3 ENSG00000117525

I also tried using org.Hs.eg.db package but get a different problem where some of my genes return NA.
My code:
  mapIds(x = org.Hs.eg.db,
     keys = targets,
     "ENSEMBL",
     "SYMBOL",
     fuzzy = TRUE,
     multiVals = first)

And this works on some but not all the genes for some reason.
    $CYR61
[1] NA

$CTGF
[1] NA

$AMOTL2
[1] "ENSG00000114019"

$ANKRD1
[1] "ENSG00000148677"

$IGFBP3
[1] "ENSG00000146674"

$F3
[1] "ENSG00000117525"

If anyone can solve either of these problems that would be great, thanks in advance!

Have you considered just downloading a list of gene IDs and the associated gene names from Ensembl Biomart?

Yeah but I wanted to be able to change out the genesets quickly and this just added another annoying step

There's an R package to allow fetching those programmatically from biomart.

You can pull a full list of Ensembl gene IDs and MGI symbols from biomaRt in R. From there it should be easy enough to filter for the genes you want:
library(biomaRt)
hsapiens_genes <- getBM(attributes = c("ensembl_gene_id",
                                       "hgnc_symbol"),
                        mart = useMart("ensembl", dataset = "hsapiens_gene_ensembl"))

You do not find some genes in your example because you were not using the official gene name. For example CTGF is offocially called CCN2 and with this name you will find the appropriate Ensembl gene ID in the biomaRt table. This is why I personally prefer to use gene name that are like ENSG00000118523|CCN2 in my work. That way you always have a unique identifier (the Ensembl ID) that is unambiguous and you do not run into trouble that for one gene many more or less official gene names exist.


text:
DESeq2 compare within a condition
This might be a really stupid question, but I can't figure out how to do this. I've read through the DESeq2 vignette and manual pages but couldn't find an answer.
I have a bunch of samples split up into different conditions (eg. celltypes and disease state). I would like to make a comparison between the two possible values within one condition, but only for the samples with a specific value of the other condition.
For example, with the following samplesheet:
patient | phenotype | type
--------+-----------+-----
1       | healthy   | A
1       | healthy   | B
1       | sick      | A
1       | sick      | B
2       | healthy   | A
2       | sick      | A
2       | sick      | B

I would like to compare "healthy" vs "sick", but only for type "A".
Currently, I have the following code, but this will also include the type "B" sample:
dds <- DESeqDataSetFromMatrix(countData=counts, colData=design, design = ~ patient + phenotype + type)
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
dds <- DESeq(dds)
res <- results(dds, contrast=c("phenotype", "healthy", "sick"))

Any idea how to accomplish this?

Welcome to the world of all possible pairwise comparisons in DESeq2 !
So the easiest way is to create another group in your data frame, first I simulate some sensible counts:
counts = counts(makeExampleDESeqDataSet(m=48))
design = expand.grid(id=1:3,
phenotype=rep(c("sick","healthy"),each=2),
type = rep(c("A","B"),2))

And you make a group that is a combination of phenotype and type:
design$group = paste(design$phenotype,design$type,sep="_")

Run DESeq2, including id:
dds <- DESeqDataSetFromMatrix(countData=counts, colData=design,
design = ~ id + group)
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
dds <- DESeq(dds)
res <- results(dds, contrast=c("group", "healthy_A", "sick_A"))

So this compares only healthy_A with sick_A.. hopefully this is what you need... The other option (normal statistical) is to fit a model with interaction, but I think this is easier to explain and interpret.
Another question you might have is why not do two separate analysis for A and B, I think it's sometimes better if you model the variance with all available data, given the limited number of replicates.


text:
Is it possible to go from a digital sequence to an actual nucleotide?
I'm trying to learn about what the state of technology is at presently.   It seems like we clearly can go from nucleotide to digitally stored sequence, but can we transcribe something from the digital sequence?

Are you asking about synthetic biology? Can you give a bit more context? The answer is it depends, e.g. one can trivially synthesise short primers but not a whole human genome

@Chris_Rands, you mean an amount of material required for an insertion is doable, but you can't synthesize a complete gene?  I'm trying to understand the technological state of bio-engineering, and trying to understand the workflow.

Oligonucleotide synthesis is a service offered by many providers such as IDT.

What is interesting is how phosphoramidite synthesis differs from how polymerases do it. It is a solid state synthesis approach that expands from the 3' to the 5' in cycles of reaction and 2'-deprotection to allow a single nucleotide to be added at a time.

To add to @Maximilian Press' answer.
There is quite a lot of research on how to store arbitrary digital data in DNA. Nick Goldman describes his research on the topic in his TEDx talk. The last progress in the field I have noticed was DNA fountain, a tranformation of bites in DNA with maximal efficiency of data storage in DNA (information / nt) and reliable recovery when sequenced afterwards.


text:
CCP4 file to a Python 3 numpy array or similar workaround
I would like to merge together several ccp4 formatted density maps (and do a few minor things).
So ideally I would like to open the ccp4/mrc files as numpy arrays in Python 3 and save the array as a ccp4 file. CCTBX can open ccp4 maps as numpy arrays, but I really would like to avoid using Python 2.7 for various reasons.
This is to be used for constraining energy minimisation in PyRosetta using the pyrosetta.rosetta.protocols.electron_density.SetupForDensityScoringMover mover/sampler. But its pyrosetta.rosetta.core.scoring.electron_density.getDensityMap data loader only likes ccp4/mrc files (density maps) and not mtz files (structure factor maps).
From what I can tell there is no way to fudge a ccp4 map from a mtz map, correct? The fact that only ePDB provides ccp4 maps makes me want to use mtz files.
CCP4's Blend can merge datasets but uses only mtz files...

You can use Gemmi.
import gemmi
ccp4_map = gemmi.read_ccp4_map('my.map')
ccp4_map.setup()  # optional
np_array = numpy.array(ccp4_map.grid, copy=False)
# ...
ccp4_map.write_ccp4_map('out.ccp4')

The setup() call expands data to the whole unit cell.
The numpy array above accesses the data through Python's Buffer Protocol (it does not copy the data). So changing values in numpy array changes gemmi.Ccp4Map.
You can also use gemmi to read mtz file and transform map coefficients to a real-space map, see https://gemmi.readthedocs.io/en/latest/hkl.html.

Two birds in one stone! Thanks ever so much!


text:
phmm blast hmmer
i have created a phmm in hmmer from a multiple alignment focused on E. coli produced in tcoffee from sequences taken from blast similar to my query sequence.

I would like to ask how i can compare the results from HMMER for a pHMM against UNIPROT (i used the phmm that i created and i searched in hmmer search against uniprot) and the results from BLAST.
I would like to ask if i can search a pHMM in BLAST somehow, or if blast accepts only fasta one sequence.
Why e-values of hmmer and blast are not comparable??Why e-value in hmmer is never 0?
Thank you so much!!

It would be useful to add the biological system you are working on. Multi-part questions don't really work, so the second part of your question is a separate question in my opinion. 1. Practicality of HMMER->UNIPROT. 2. Theories of e-values.  Blast only accepts fasta, but conversion is easy, so programmatically its doable.

thanks!yes i should create a second post. it is about E.coli. and i am trying to compare the results taken from hmmer with the results taken from blast!thanks!


text:
how can i search in p-blast against UNIPROT(not only swissprot)
i need to search for one protein against UNIPROT in BLAST for a selected organism.
How can I do it? (SWSISSPROT+TREMBL)
THANKS!!

I don't know of any online tool enabling a search against TrEMBL. You probably will have to set up your own local BLAST environment, roughly following these steps:

download blast+ - but there might be an easier way to install the software, depending on your operating system (e.g. on Ubuntu you can just sudo apt install ncbi-blast+)
download target databases
build local database

I hope this helps.

thank you a lot for your answer!unfortunately my programming skills are not so strong!

UniProt provides a BLAST service (https://www.uniprot.org/blast) which you can run against several UniProt databases, or you could use the Expasy BLAST too, at https://web.expasy.org/blast.