cannin/biostars_qa huggingface dataset:
Question: So I was struggling with this: Creating a dendrogram with a large dataset (20,000 by 20,000 gene-gene correlation matrix): https://www.biostars.org/p/477332

Now I think I found a good solution, taking a "meaningful" sample of the dataset, and then create a dendrogram-heatmap of the gene-gene correlation matrix generated from the sample.

I have got this far:

cluster3.seurat.obj <- CreateSeuratObject(counts = cluster3.raw.data, project = "cluster3", min.cells = 3, min.features = 200)

cluster3.seurat.obj <- NormalizeData(cluster3.seurat.obj, normalization.method = "LogNormalize", scale.factor = 10000)
cluster3.seurat.obj <- FindVariableFeatures(cluster3.seurat.obj, selection.method = "vst", nfeatures = 2000)

Now I am wondering, how do I extract a data frame or matrix of this Seurat object with the built in function or would I have to do it in a "homemade"-R-way?

I'm hoping it's something as simple as doing this:

cluster3.cells.variable.features <- as.matrix(GetAssayData(cluster3.seurat.obj, slot = "data")[, WhichCells(cluster3.seurat.obj)][, FetchData(cluster3.seurat.obj, var.features)])

I was playing around with it, but couldn't get it...

Any help would be appreciated.

Very Respectfully,
Pratik

Answer: You just want a matrix of counts of the variable features?

var_genes <- VariableFeatures(cluster3.seurat.obj)
seurat_df <- GetAssayData(cluster3.seurat.obj)[var_genes,]
Question: Hello,

I have a problem when using Music. Only missense, nonsense, splicing mutations and indels are included in my data, there are no synonymous mutations, but the tumor do have synonymous mutations.

Is it OK to find SMGs with Music ? Is it OK to define background mutation without silent mutations?

Thank you~

Answer: Read the arguments of `music bmr calc-bmr`, where it defaults to excluding Silent mutations. The documentation is [online here](http://gmt.genome.wustl.edu/packages/genome-music/genome-music-bmr-calc-bmr.html). If you don't want MuSiC to skip Silent mutations, use the argument `--noskip-silent`.

But it is preferable to `--skip-silent` when finding SMGs with somatic mutations, because silent mutations are more likely to be passenger events with lower selective pressure during cancer evolution (though there are exceptions, like when Silent mutations alter promoter motifs). You can try running it with `--skip-silent` and `--noskip-silent`, and compare results.
Question: I am currently working with a program that generates aritificial FASTQ files when given a reference genome called Artificial FASTQ Generator. [Here][1] is a link to the description of the program and [here][2] is the manual, it says that the program generates paired-end reads (and it does generate two FASTQ files). After aligning the artificial reads (the two FASTQ files generated by Artificial FASTQ Generator) to the reference genome as paired-end reads using Bowtie 2 I got the following result:

```
892589 (100.00%) were paired; of these:
892585 (100.00%) aligned concordantly 0 times
4 (0.00%) aligned concordantly exactly 1 time
0 (0.00%) aligned concordantly >1 times
----
892585 pairs aligned concordantly 0 times; of these:
870486 (97.52%) aligned discordantly 1 time
----
22099 pairs aligned 0 times concordantly or discordantly; of these:
44198 mates make up the pairs; of these:
7154 (16.19%) aligned 0 times
15716 (35.56%) aligned exactly 1 time
21328 (48.26%) aligned >1 times
99.60% overall alignment rate
```

After aligning both FASTQ files as single end reads, I got the following:

```
1785178 (100.00%) were unpaired; of these:
7154 (0.40%) aligned 0 times
1755694 (98.35%) aligned exactly 1 time
22330 (1.25%) aligned >1 times
99.60% overall alignment rate
```

What I do not understand is why these reads are aligning as single-end reads and not as paired-end reads as expected? Is anybody familiar with both programs that can help explain this?

EDIT:

Here are the outputs:

- [1st FASTQ output of AFG][3]
- [2nd FASTQ output of AFG][4]
- [SAM output of Bowtie 2][5] (single-end alignment)

EDIT 2:

Here is the output of the paired-end alignment.

[SAM output of Bowtie 2][6] (paired-end alignment)

[1]: http://sourceforge.net/projects/artfastqgen/
[2]: http://sourceforge.net/p/artfastqgen/wiki/Home/
[3]: https://drive.google.com/file/d/0B3xfexs8wgRTS2tiYkROanJYdnM/edit?usp=sharing
[4]: https://drive.google.com/file/d/0B3xfexs8wgRTRWdxeV9iZkhoOGc/edit?usp=sharing
[5]: https://drive.google.com/file/d/0B3xfexs8wgRTOVV2Rl9KSUxSc3c/edit?usp=sharing
[6]: https://drive.google.com/file/d/0B3xfexs8wgRTR1NDM3A3TE5JX2c/edit?usp=sharing

Answer: Maybe you should post an example of the output, a few read pairs should be sufficient.

There could be many reasons, e.g.:

- Bug or wrong use of the simulator: mate pairs not on complementary strands
- too large insert size
- Bug: annotated pairs do not pair
- too low simulated quality scores (single-end do align though, so not in this case)
- too high simulated error rate (dto.)

Maybe you should use a simulator that is better tested, see:

https://www.biostars.org/p/2194/
Question: All proteins on UniProt have a unique accession number. Ex "O15169" is the accession for human Axin 1.

Other RDF stores referring to proteins on UniProt use this accession (eg [Pathway Commons reference][2])

[This document][1] describes the RDF schema for UniProt.

Where is the UniProt accession in this RDF schema?

[1]: https://www.uniprot.org/core/
[2]: http://rdf.pathwaycommons.org/describe/?url=http://pathwaycommons.org/pc12/UnificationXref_uniprot_knowledgebase_O15169

Answer: In the UniProt RDF model, the accession is only in the IRI of the form `http://purl.uniprot.org/uniprot/${ACCESSION}`.

To go from an accession string in pathway commons to a IRI one uses a SPARQL snippet like:

VALUES ?acc { "P05067" }
BIND(IRI(CONCAT("http://purl.uniprot.org/uniprot", ?acc)) AS ?entry)

There are two reasons that we don't have the primary accession as a string in our RDF or SPARQL endpoint.

1. Avoiding false joins, an UniProt accession. Might also be used to identify something completely else, without the IRI part false joins can lead to wrong results.
2. Adding a string for each identifier adds hundreds of millions of extra triples and strings in the database which will negatively impact performance and storage.
Question: I have a list of accession numbers, how can I search for there function in Pfam db http://pfam.xfam.org/ as a bulk?

For example:

If I have this accession **PF07714.12** and I searched manually I will get;

> A tyrosine kinase is an [enzyme][1] that can transfer a [phosphate][2] group from [ATP][3] to a [protein][4] in a cell. It functions as an "on" or "off" switch in many cellular functions. Tyrosine [kinases][5] are a subclass of [protein kinase][6].
>
> The phosphate group is attached to the [amino acid][7] [tyrosine][8] on the protein. Tyrosine kinases are a subgroup of the larger class of [protein kinases][9] that attach phosphate groups to other amino acids ([serine and threonine][10]). [Phosphorylation][11] of proteins by kinases is an important mechanism in communicating signals within a cell ([signal transduction][12]) and regulating cellular activity, such as cell division.
>
> Protein kinases can become mutated, stuck in the "on" position, and cause unregulated growth of the cell, which is a necessary step for the development of cancer. Therefore, kinase inhibitors, such as[imatinib][13], are often effective cancer treatments.

How can I do it for a list of accessions?

[1]: http://en.wikipedia.org/wiki/Enzyme
[2]: http://en.wikipedia.org/wiki/Phosphate
[3]: http://en.wikipedia.org/wiki/Adenosine_triphosphate
[4]: http://en.wikipedia.org/wiki/Protein
[5]: http://en.wikipedia.org/wiki/Kinases
[6]: http://en.wikipedia.org/wiki/Protein_kinase
[7]: http://en.wikipedia.org/wiki/Amino_acid
[8]: http://en.wikipedia.org/wiki/Tyrosine
[9]: http://en.wikipedia.org/wiki/Protein_kinase
[10]: http://en.wikipedia.org/wiki/Serine/threonine-specific_protein_kinase
[11]: http://en.wikipedia.org/wiki/Phosphorylation
[12]: http://en.wikipedia.org/wiki/Signal_transduction
[13]: http://en.wikipedia.org/wiki/Imatinib

Answer: Unlike what you have shown here, if **one sentence description** is enough to describe about the family you can use `hmmfetch` which first extract profiles and from that file you can separate ACC & DESC lines easily, ACC - accession id, DESC - description about family. [Here][1] you can download all pfam family HMM profiles in one file, then use it as a source to extract families of your interest.

[1]: ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release
Question: Hello,

I want to extract some specific strings in csv file, using either excel, R or Python.

for example as below: I want to find string from column A in column B and return in column C with 5 amino acid before and after **N**; thanks!!

```
A B C
INETTDFR MHRFLLMLLFPFSDNRPMMFFRSFIVFFFLIFFASNVSSRKQTYVIHT IVGKINETTDF
VTTSTKHIVTSLFNSLQTENINDDDFSLPEIHYIYENAMSGFSATLTDDQLDT
VKNTKGFISAYPDELLSLHTTYSHEFLGLEFGIGLWNETSLSSDVIIGLVDTG
ISPEHVSFRDTHMTPVPSRWRGSCDEGTNFSSSECNKKIIGASAFYKGYE
SIVGKINETTDFRSTRDAQGHGTHTASTAAGDIVPKANYFGQAKGLASGM
RFTSRIAAYKACWALGCASTDVIAAIDRAILDGVDVISLSLGGSSRPFYVDP
IAIAGFGAMQKNIFVSCSAGNSGPTASTVSNGAPWLMTVAASYTDRTFPAIV
RIGNRKSLVGSSLYKGKSLKNLPLAFNRTAGEESGAVFCIRDSLKRELVEGK
IVICLRGASGRTAKGEEVKRSGGAAMLLVSTEAEGEELLADPHVLPAVSLGF
SDGKTLLNYLAGAANATASVRFRGTAYGATAPMVAAFSSRGPSVAGPEIAKP
DIAAPGLNILAGWSPFSSPSLLRSDPRRVQFNIISGTSMACPHISGIAALIKSV
HGDWSPAMIKSAIMTTARITDNRNRPIGDRGAAGAESAATAFAFGAGNVDPT
RAVDPGLVYDTSTVDYLNYLCSLNYTSERILLFSGTNYTCASNAVVLSPGDLN
YPSFAVNLVNGANLKTVRYKRTVTNVGSPTCEYMVHVEEPKGVKVRVEPKVL
KFQKARERLSYTVTYDAEASRNSSSSSFGVLVWICDKYNVRSPIAVTWE
```

Answer: Using Excel*:

=MID(B1,FIND(A1,B1)+FIND("N",A1)-6,5) & "N" & MID(B1,FIND(A1,B1)+FIND("N",A1),5)

*Don't use Excel :)
Question: Hey everyone,

I'm inexperienced with statistics, and want to perform a regression between two diseases. I would really appreciate some clarification if my understanding of when to include an interaction term is correct.

Let's call the diseases disease Y and disease X. I know that disease Y is age-dependent, i.e. it's more likely to be encountered in older individuals and it becomes progressively worse with age. So I think the regression should definitely include age as a covariate:

Y ~ age + X

I'm not sure whether disease X is also age-dependent, but it might be. My planned approach was to look into the data, check if an independent student's t-test detects a significant difference in the age distributions between people with / without disease X. If yes, I would correct the regression formula to include an interaction term:

Y ~ age + age:X + X

Would this approach be correct?

Additionally, would it matter if my variable Y represents case/control status or disease severity (i.e. logistic vs linear regression)?

Answer: It's best practice to first check if your variables are correlated. If they are, you should either drop one or combine them into one variable.
In R:

cor.test(your_data$age, your_data$X)

I would drop one of the variables if r >= 0.5, although others may use a different cutoff. If they are correlated, I would keep the variable with the lowest p-value. Alternatively, you could combine age and X into one variable by adding them or taking their average.
To find p-values:

model = lm(Y ~ age + X, data = your_data)
summary(model)

If age and X are not correlated, then you can see if there is an interaction.

int.model = lm(Y ~ age + X + age:X, data = your_data)
summary(int.model)

If the interaction term has a significant p-value, then you'll want to include it in your model. If not, then you'll want to drop it.
You can use either linear or logistic regression. For logistic regression, you would use the following:

logit.model = glm(Y ~ age + X + age:X, data = your_data, family = binomial)
summary(logit.model)
Question: Do you know any public scientific SQL server ?

for example, I would cite:

UCSC http://genome.ucsc.edu/FAQ/FAQdownloads#download29

ENSEMBL http://uswest.ensembl.org/info/data/mysql.html

GO http://www.geneontology.org/GO.database.shtml#mirrors

(I'll give a +1 to each correct answer)

Answer: This is quite an important one for people doing mouse work, though it is important to note that JAX offer Mart and Batch-Query functionality through the web site as well which may well suit many peoples needs.

- Direct SQL Access to MGI
- http://www.informatics.jax.org/software.shtml#sql

Note that this is a public 'free' service, but that you do need to contact user support to get your login and password. They are also happy to provide some custom SQL scripts to get you started.
Question: Heys,
I'm working with bedtools trying to find out the coverage of several target regions.
I have my bam files and my target regions in a .bed file and I'm running this:

bedtools bamtobed -i WA01.bam | bedtools coverage -a - -b mit.bed > exons.mit.coverage

When I open the exons.mit.coverage file, in the coverage section I just have 0s in the areas that are not matching my target genes and 1s in the areas of my target genes, as I have **coverage per read**.

My question: How can I obtain coverage per gene and not coverage per read?
I imagine I can filter this out, selecting the reads with a good quality score per gene, but I think this should be already implemented anywhere else that I don't know

Thanks for your help! :)

Answer: Switch `-a` and `-b`.

bedtools bamtobed -i WA01.bam | bedtools coverage -b - -a mit.bed > exons.mit.coverage
Question: According to http://www.ebi.ac.uk/~zerbino/velvet/, velvet can be multithreaded. But how to do it? When I do

```
./velveth

Usage:
./velveth directory hash_length {[-file_format][-read_type][-separate|-interleaved] filename1 [filename2 ...]} {...} [options]
```

My question is how to run velveth multithreaded? Where is the multi-thread flag (usually `-t` or `-p`)? I'm running 1.2.10.

Answer: For threading to work, Velvet has to be installed with make 'OPENMP=1'. Then you can control threading via export OMP_NUM_THREADS=7. See Velvet Manual Section 2.3.6.
Question: As per the title: why?

Does the preliminary PCA step merely maximise the chances of identifying the clusters that tSNE / UMAP later identify?

When doing this, in tutorials, it seems that people blindly choose a number of PCs as input to tSNE / UMAP, without checking the amount of variation explained by the PCs. For example, 20 PCs may only account for 30% overall variation.

Assuming one follows a standard procedure for, e.g., scRNA-seq, is this one step too many. A pipeline could be:

1. normalisation of raw counts
2. transformation of normalised counts
3. PCA-transformation (may include an additional pre-scaling step)
4. tSNE/UMAP

Kevin

Answer: At least from a clustering perspective, I'd probably try it both ways to be on the safe side, i.e. with PCA to get the top N PCs and without. I'm a bit skeptical of reducing to N PCs for clustering because there is inevitable information loss. The same will apply for t-SNE, UMAP, etc. I'd prefer to use the most variable genes instead.

Although I think it is less of a problem when just trying to visualise the data, rather than define some new cell types using cluster analysis where we might want to be a bit more cautious.
Question: Hello,

I have multiple fasta sequences that are like this:

>2p__scaffold_2__5799__6580__-__778568__0.00__0.00
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAACTGAGAGTCGGCGATGAGCAGCCGGA
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAA
>2p__scaffold_2__5799__6580__+__778569__0.00__0.00
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAACTGAGAGTCGGCGATGAGC
>1p__scaffold_2__11235__11438__-__830827__0.00__0.00
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAACTGAGAGTCGGCGATGAGCAGCCGGA
GCTTCAATCCAGGGGATCGAGGAGATCCAAAGCAGCAGAAGCGGCTCGACGATGGTGAGG
ATTCGGGATCGGATTCAGCGCTCGTCGGGACTGG
>1p__scaffold_2__33129__34129__+__811706__0.00__0.00
GCTGGCGACGGATCTA

And I want to keep just the "> + ID" (numbers after `__+/-__` and before __0.00_0.00)

So I expect an output like this:

>778568
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAACTGAGAGTCGGCGATGAGCAGCCGGA
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAA
>778569
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAACTGAGAGTCGGCGATGAGC

I searched for it and tried this:

sed 's@.*__-__@@' input.fa > output.fa

That removed `__-__` and everything before it, including the ">" that I wanted to keep.

I also tried this to remove everything between ">" and `__-__`

sed -e 's/\>//' -e 's/\__-__.*//' input.fa > output.fa

But this removed everything after `__-__`

And this, that removed __0.00_0.00

sed 's/__0.00.*$//' input.fa > output.fa

Thank you for your help.

Answer: This should work if your headers follow the pattern that you specified:

sed 's/[>_]\+/_/g' yourfile.fasta | cut -f 8 -d _ | sed 's/^\([0-9]\)/>\1/'
Question: Hello,

I want to extract some specific strings in csv file, using either excel, R or Python.

for example as below: I want to find string from column A in column B and return in column C with 5 amino acid before and after **N**; thanks!!

```
A B C
INETTDFR MHRFLLMLLFPFSDNRPMMFFRSFIVFFFLIFFASNVSSRKQTYVIHT IVGKINETTDF
VTTSTKHIVTSLFNSLQTENINDDDFSLPEIHYIYENAMSGFSATLTDDQLDT
VKNTKGFISAYPDELLSLHTTYSHEFLGLEFGIGLWNETSLSSDVIIGLVDTG
ISPEHVSFRDTHMTPVPSRWRGSCDEGTNFSSSECNKKIIGASAFYKGYE
SIVGKINETTDFRSTRDAQGHGTHTASTAAGDIVPKANYFGQAKGLASGM
RFTSRIAAYKACWALGCASTDVIAAIDRAILDGVDVISLSLGGSSRPFYVDP
IAIAGFGAMQKNIFVSCSAGNSGPTASTVSNGAPWLMTVAASYTDRTFPAIV
RIGNRKSLVGSSLYKGKSLKNLPLAFNRTAGEESGAVFCIRDSLKRELVEGK
IVICLRGASGRTAKGEEVKRSGGAAMLLVSTEAEGEELLADPHVLPAVSLGF
SDGKTLLNYLAGAANATASVRFRGTAYGATAPMVAAFSSRGPSVAGPEIAKP
DIAAPGLNILAGWSPFSSPSLLRSDPRRVQFNIISGTSMACPHISGIAALIKSV
HGDWSPAMIKSAIMTTARITDNRNRPIGDRGAAGAESAATAFAFGAGNVDPT
RAVDPGLVYDTSTVDYLNYLCSLNYTSERILLFSGTNYTCASNAVVLSPGDLN
YPSFAVNLVNGANLKTVRYKRTVTNVGSPTCEYMVHVEEPKGVKVRVEPKVL
KFQKARERLSYTVTYDAEASRNSSSSSFGVLVWICDKYNVRSPIAVTWE
```

Answer: Something like this should get you started. I haven't tested this, though, and you may need to change `delimiter` and `quotechar` parameters based on what your input and output look like:

```
#!/usr/bin/env python

import csv
import sys

padding = 5
assert padding >= 0, 'padding is less than zero, which does not make sense'

residues_of_interest = 'N'
roi_length = len(residues_of_interest)

csv_writer = csv.writer(sys.stdout, delimiter=',', quotechar='\"', quoting=csv.QUOTE_ALL)
for row in csv.reader(iter(sys.stdin.readline, ''), delimiter=',', quotechar='\"'):
A = row[0]
B = row[1]
pos_roi_in_A = A.find(residues_of_interest)
if pos_roi_in_A == -1:
continue # skip over lines where residues-of-interest string is not in A
pos = B.find(A) + pos_roi_in_A # anchor 'pos' value on the start position of ROI
lower_bound = padding
upper_bound = len(B) - padding
#
# consider four cases of a bounded string search result
#
if pos == -1:
C = None
elif pos < lower_bound:
C = B[0:(pos + roi_length + padding)]
elif pos > upper_bound:
C = B[(pos - padding):upper_bound]
else:
C = B[(pos - padding):(pos + roi_length + padding)]
#
# write A, B and C to standard output
#
csv_writer.writerow([A, B, C])
```

You might run it like so:

$ padded_substring_finder.py < in.csv > out.csv
Question: hi friends,

i RNA-seq protocol i read...

tophat -p 8 -G genes.gtf -o C1_R1_thout genome C1_R1_1.fq C1_R1_2.fq

by which i have accepted_hits.bam

cufflinks -p 8 -o C1_R1_clout C1_R1_thout/accepted_hits.bam

then saying that :

> Create a file called assemblies.txt that lists the assembly file for each sample. The file should contain the following lines:
./C1_R1_clout/transcripts.gtf
./C2_R2_clout/transcripts.gtf
./C1_R2_clout/transcripts.gtf
./C2_R1_clout/transcripts.gtf
./C1_R3_clout/transcripts.gtf
./C2_R3_clout/transcripts.gtf

then this command

cuffmerge -g genes.gtf -s genome.fa -p 8 assemblies.txt

what is the `assemblies.txt` please??? how i can create such a file by *cufflinks output* that as i saw is a file named `QC_filtered.fastq`

thanks for your simpathy

Answer: The output of cufflinks is transcripts.gtf . If you have many samples, cufflinks creates transcripts.gtf ( i.e assembles transcripts) file for each samples and these needs to be merged by cuffmerge,
Question: Hi all,

I need to assign GO terms to DE genes in R. `ClusterProfiler` seems to be a popular package for that.
But for some reason it either fails to detect GO terms or its GO terms do not coincide with other tools.

library(clusterProfiler)
library(org.Hs.eg.db)

###mtDNA genes
genes<-c("ENSG00000198886","ENSG00000198888","ENSG00000198763","ENSG00000198840","ENSG00000212907","ENSG00000198786","ENSG00000198695","ENSG00000198727","ENSG00000198804","ENSG00000198712","ENSG00000198938","ENSG00000198899","ENSG00000228253")

ego<-enrichGO(genes, OrgDb = org.Hs.eg.db, keytype ="ENSEMBL",ont = "CC",pAdjustMethod = "fdr",pvalueCutoff = 0.1, qvalueCutoff = 0.5)

dotplot(ego, showCategory=30)

The same genes with `Panther` as expected show mtDNA related GO terms. Other genes from my dataset get enrichments, so apparently the problem is not in gene IDs.
What am I doing wrong?

Cheers,

Answer: With the hint from @Guangchuang in the post above, I resolved the problem by updating `clusteProfiler` to the latest version.

For that I had to update R to v. 3.5 for my Ubuntu 16.04 https://www.r-bloggers.com/updating-r-on-ubuntu/

Then install the latest bioconductor version https://www.bioconductor.org/install/

And resolve the igragh fortrat4 dependancy issue https://ashokragavendran.wordpress.com/2017/10/24/error-installing-rigraph-unable-to-load-shared-object-igraph-so-libgfortran-so-4-cannot-open-shared-object-file-no-such-file-or-directory/

Hope this will save several hours for people with the same problem
Question: Hi, I'm pretty new to the scRNA-seq world and while working on my own sets of data, I'm starting to wonder when should the batch correction algorithm be used appropriately.

Let's say we have a Day0, a Day1, a Day2 and a Day3 scRNA-seq sample.

To elaborate, starting from Day0, assume we treated certain chemical, and sampled it on a daily basis during the course of experiment.

**Would it be ok or reasonable to apply batch correction algorithm (e.g. CCA) to this aggregation of samples?
I mean, is CCA algorithm designed for this kind of experiment design?**

From the experiment from [Kang et al., 2017][1] which is comprised of PBMC, splitted into a control group and a stimulated group treated with interferon beta, [they state that][2] "the repsonse to interferon caused cell type specific gene expression changes that makes a joint analysis of all the data difficult with cells clustering both by stimulation condition and by cell type". But is it reasonable?

My understanding is that if you are to use batch correction you should have biological or technical batches from the "same condition". So if you have replicate samples with the same condition and when somehow they are separated from each other for technical reason, it's appropriate to use batch correction.

Going back to the supposed experiment I stated above, I think (maybe i'm wrong and i am most of time) it's not reasonable to apply batch correction to this Day0-4 experiment.

Can someone give me some clear explanation to the use of batch correction?

Thank you.
Ryan

[1]: https://www.nature.com/articles/nbt.4042
[2]: https://satijalab.org/seurat/v3.0/immune_alignment.html

Answer: > if you are to use batch correction you should have biological or technical batches from the "same condition"

From the [integration vignette][1]: "These methods aim to identify shared cell states that are present across different datasets, even if they were collected from different individuals, **experimental conditions**, technologies, or even species"

Thus, the "official" answer is that different conditions are fine.

Really, it depends on the questions you want to ask and on the data that you have. For example, if all your time points segregate and form distinct clusters, it's going to be hard to present any kind of coherent analysis.

[1]: https://satijalab.org/seurat/v3.0/integration.html