cannin/biostars_qa huggingface dataset:
Question: Hi,

So I have two datasets from two different yet related cell lines: pre and post relapse cancer cell line from same patient.

I have performed single cell sequencing for both of them with the hypothesis that I will be able to find rarely expressed cells in either that are key to becoming relapse.

To do that I am following two protocols:

1) I am following the Seurat tutorial of integrating simulated and normal pbmc

https://satijalab.org/seurat/v3.1/immune_alignment.html

2) I am following this second tutorial where there is no integration steps involved.

https://davetang.org/muse/2018/01/24/merging-two-10x-single-cell-datasets/

However both these approaches give two different results! Using 1) I am getting more cells from two different lines to be similar with very few cells from each to be different. Using 2) I am getting opposite, yet expected, result that most of the cells are different with only few common cells.

Which one is the correct way of analyzing? Are two approaches giving me two different results become they are doing different things? 1) is finding common genes and 2) finding distinct genes?

Answer: These are very different approaches. The first tries to account for technical variation between the two sets, mapping similar cells between the two to each other. The second literally just merges the columns and rows of the two sets together into one object - it's not doing any special normalization. It's just a straight merging of data.

It's tough for us to say which is more appropriate - the first may help you to identify cell populations that truly differentiate the two, but it could be blowing away real differences due to how Seurat's integration works. It will force populations that aren't similar together if there aren't many overlapping cell types between the two samples. The second may be revealing significant technical variation or batch effects, or it could just be that your cell lines are quite different from each other. You are in the best position to determine if this is the case or not - we know nothing about your samples.

This is where the true difficulty of RNA-seq analysis lies - nobody is really going to be able to tell which is truly correct.

You might try other integration methods if you feel you have batch effects or technical variation that needs to be addressed. I've found the `SeuratWrapper` around [fastMNN](https://htmlpreview.github.io/?https://github.com/satijalab/seurat.wrappers/blob/master/docs/fast_mnn.html) to be quite good, personally, as it handles cases where samples don't have much overlap in terms of cell types much more appropriately.
Question: The "snp" command works with only 1 snp, but if I type 2 SNPs, it does not, as from output below. VCF tools manual says snp but how do I write the string of several SNPs? rs236330rs236330 or rs236330-rs236330?

I tried the command "snps" to include a list of SNPs from a file but it doesn't recognize the command

```
--vcf ALL.chr1.phase3_shapeit2_mvncall_integrated_v5.20130502.genotypes.vcf
--chr 1
--freq
--out chr1_analysis
--snp rs236330-rs236322

After filtering, kept 2504 out of 2504 Individuals
Outputting Frequency Statistics...
After filtering, kept 0 out of a possible 6468094 Sites
No data left for analysis!
```

Answer: I believe you can simply use `--snp` multiple times, one for each snp. Alternatively, you can read from a file using `--snps`.
Question: Using LIMMA's removeBatchEffect and PVCA to check for batch influences before and after removal I just noticed that if I don't use log2 values the batch effect removal does not work in the sense that PVCA shows retained batch effects after removal. With log2 values all runs fine.

Is this a general feature of removeBatchEffect?

Its a bit strange as later on, eg for VOOM, non-log2 values are required.
(though to be fair not for the corrected data as those won't go into the diff expression analysis).

Thanks for your help!
Sebastian

Answer: There are multiple tools live SVA, COMBAT, edgeR function, limma function, RUVSEq, etc.

[PVCA][1] is doing PCA + Variance component analysis to extract factors that influence your variability using linear mixed model
Steps :

1. So normalized data + PVCA = confounders of batch effects (Assessing and identifying confounders)
2. Perform correction or adjustments of the confounders extracted from 1 using any standard batcheffect adjustment methods and perform again
PVCA to visualize on normalized-adjusted- log transformed data to see if confounders are assessed and what you see as primary gnee variability is due to biological phenotypic variation.( Viewing post-batch adjustments on all the genes expressed in all samples)
3. However, for DE analysis it is the model matrix and the model
effects/covariates that will be pulled out from 1(batch effect confounders) . These effects
should be modeled around your count data for any linear model
fitting. Any results of DE should be then viewed in with log transformed/corrected (from 2) data for visualizations like heatmaps, expression box plots, etc. (using confounders as covariates in model design for linear fit with limma on counts data to perform DEA)

My only two cents are, batch effect removal is not the key, one needs to adjust for it rather not deduct it. You are trying to understand what are the confounders in your data and how they mess around. You do it from counts data and any normalization that entails should be used while performing the effect analysis. What you see for plots later can be viewed via log transformation. However, for downstream differential analysis you don't use log2 transformed batch corrected data. One uses the counts data, pulls out the effect information, adds it either as covariates in model design for differential expression or use it for adjustments. I suggest to you take a look at the below links to understand how it is done, what is the underlying statistics and variance associated , and that you do not make any over-fitting.

[Removing continuous covariate effects in limma analysis][2]

Good luck!

Edit 1: I just realized there was a bit wrong flow I suggested earlier. So corrected it.

[1]: https://www.niehs.nih.gov/research/resources/software/biostatistics/pvca/index.cfm
[2]: https://support.bioconductor.org/p/66251/
Question: More as an educational than a research tool, I'd like to give a presentation on human genetics using a human genome browser that can start at the single-base level, and zoom (smoothly if possible) out to the entire chromosome. The catch is that I won't have internet access when giving the presentation. I'm happy to download reference sequences, annotations, etc. The main things I'd like to display are (at the fine scale) intron/exon positions, and at the wider scale, gene positions, alu/SINE positions, and (if possible) simple gene repeats such as human MW/LW opsins. I like http://chromozoom.org, but it's hard to get working offline. I find that http://www.biodalliance.org doesn't have a terribly nice chromosome-wide view. Are there others that people could suggest?

Answer: [Integrated Genome Browser [IGB]][1] may be what you are looking for. We used to use this for research as well but at the time it did not scale well for large datasets (I have not performance tested it recently though)

[Integrative Genomics Viewer [IGV]][2] is what we use at the moment for research but it might not look as slick for your intended audience.

Both genome browsers can work offline with local files (e.g. bed format).

[1]: http://bioviz.org/igb/
[2]: http://www.broadinstitute.org/igv/
Question: I am a student and I have a fun problem: My friend doesn't like cats and thinks raccoon's are "way cooler." I want to build a phylogenetic tree and show him that we're more closely related to cats and therefore cats are cooler.

However, my school assignments always told me which gene or protein to use for analyses.

So how do I decide which gene/prot to compare without looking back at old assignments?

I'll be looking at house pets so I know I need something that cats, dogs, lizards, gold fish, and snakes have in common.

What should be my next steps?

Answer: Phylogenies using a single gene are not reliable. For this kind of comparison, you should probably use fast evolving genes such as mitochondrial genes. Check the literature, people have most likely already worked out which genes to use and how.
Question: Hi,
I want to perform differential expression analysis of a GSE Series **GSE138214** (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE138214) using limma in R. I am unable to extract data from this series as this series has around 14 different platforms.

Using GEOquery and limma package, I am able to differentially analyse genes in Series **GSE138260** (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE138260), since it has only one platform.

But the GSE138214 series has 14 platforms and I am not sure how to work around it.

Thanks

Answer: It looks like you are interested in differential gene expression. Note that there are only two samples, so statistics won't work for you here. Also, these are tiling arrays, not traditional gene expression arrays. Finally, the GEO submission did not include values, so `getGEO` will have to be augmented with custom code and accessing supplemental files. In short, it will be a lot of work for potentially little reward.
Question: I'm reading an interesting paper, [Malachi et al., Cell Systems 2015][1], that talks about Variant Allele Frequency (VAF). Could someone please help me understand how they calculated this value? I couldn't find it in the methods.

**edit**

For example, let's say I have a bunch of VCF files that I would like to find VAF for. How would I go about doing that? Here's an example of a VCF file:

##fileformat=VCFv4.1
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##FORMAT=
##FORMAT=
##FORMAT=
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	genotype
chr1	808922	.	G	A	249	PASS	DP=49;AB=1.0000;RD=0;AD=49;SAP=0.5510;RAP=0;DP4=0,0,27,22	GT:GQ:DP	1/1:145:49
chr1	808928	.	C	T	249	PASS	DP=51;AB=1.0000;RD=0;AD=51;SAP=0.5294;RAP=0;DP4=0,0,27,24	GT:GQ:DP	1/1:145:51
chr1	876499	.	A	G	217	PASS	DP=49;AB=1.0000;RD=0;AD=49;SAP=0.6939;RAP=0;DP4=0,0,34,15	GT:GQ:DP	1/1:126:49

[1]: http://www.cell.com/cell-systems/fulltext/S2405-4712(15)00113-1

Answer: Variant allele frequency in this case refers to the fraction of sequencing reads overlapping a genomic coordinate that support the non-reference (mutant/alternate) allele.

Typically, this information is either explicitly listed or readily extractable from VCF files. If not, given that you have the genomic coordinates of the mutation of interest and the bam file, you can run [samtools mpileup][1] command to get the alleles coming from the entire set of overlapping reads.

I use a [python script][2] to parse the output of mpileup into a human-readable allelic table.

[1]: http://samtools.sourceforge.net/mpileup.shtml
[2]: https://github.com/Niknafs/NGSTools/blob/master/baseParser.py
Question: Hello there,

I am having a lot of MiSeq data I am trying to analyze and I figured out using FastQC that I have a lot of fails in the report it generates and I wanted to pick your brain to get sense of what should be done in that case.

As you all know FastQC generates this kind of information:

```
[PASS] Basic Statistics
[FAIL] Per base sequence quality
[PASS] Per sequence quality scores
[FAIL] Per base sequence content
[FAIL] Per base GC content
[WARNING] Per sequence GC content
[PASS] Per base N content
[WARNING] Sequence Length Distribution
[FAIL] Sequence Duplication Levels
[WARNING] Overrepresented sequences
[FAIL] Kmer Content
```

The issue here is that I am analyzing targeted sequencing data, so I am expecting to have a lot of duplications, what I don't clearly see is whether or not to take the result of FastQC as correct based on the standard they are publishing on their website (how a good report should look like and how a bad one should look like), so I am expecting the GC content to go crazy with the amount of duplication because of the type of experiment (deep sequencing), now based on the information provided in the example above, do you think fastq post processing like clipping and trimming would correct the reads or is it failing in the level of the MiSeq machine already (experimental contamination?)

Rad

Answer: I would use the FastQC results directly (the plots) and comparatively (between lanes/runs of the same type of experiment) and not rely heavily on the PASS/FAIL stuff. In general, trimming can often help things, but for targeted sequencing, the aligners do a fairly good job of soft clipping adapters and low-quality bases anyway. The better metric, in my opinion, is to do alignment and see what you get. A perfectly good FastQC result can lead to very poor mapping (sequenced wrong organism, for example). Also, a pretty poor FastQC result often results in usable data. Where FastQC can be quite helpful is in determining whether sequencing the same library will be an effective way of increasing depth. If the FastQC results suggest that there are library issues, perhaps making a new library is warranted.
Question: Hello all, I am new in plink analysis

I have Map file and I know that location (start /end) of ***CD2AP ***gene on chromosome 6 has base pairs 47,477,745 to 47,627,262

How can I use physical position column in MAP file to extract SNPs of this gene

Any help will be appreciated

Answer: plink cannot process isolated `.map` files; do you have an accompanying `.ped`? (If you don't, it might be simpler to write a short script for this job.)

If you do,

plink --file [.ped/.map filename prefix] --chr 6 --from-bp 47477745 --to-bp 47627262 --write-snplist --out [output filename prefix]

should give you the list you're looking for.
Question: Many experts claim that the negative binomial distribution is better than the Poisson distribution for modeling discrete RNA-Seq data.

Is this because

The expression profiles of organisms resemble the negative binomial distribution - i.e.if I binned real gene transcription counts from an RNA-Seq experiment it would look like the plot below.

or

The sampling error of gene expression is such that the true population mean of a gene looks like the negative binomial - i.e the true mean expression level of that gene is probably more (because of the skew) than the mean expression level of a sample of reads of that gene drawn from replicates.

or are these two concepts the same thing?

Answer: Hi Jeremy.

It depends what you mean by "modeling discrete RNA-Seq data". If differential expression is your game, you need to keep in mind that we/others model the distribution of counts for a given gene across replicates; we don't model the distribution of an individual sample, as you suggest in your first point above. So for DE, we are not (explicitly) concerned with shape of the distribution of counts for the whole organism at all.

Assuming DE is your interest, the reason we choose NB over Poisson is that, in real biological applications, there is simply more variability than Poisson can explain. Poisson is a single parameter dist'n, with mean=variance. That assumption, which is really an approximation to the binomial, is suitable only for the variability associated with sampling the same DNA population (e.g. if you sequence multiple lanes of the same DNA, and assume no lane-specific effects, etc.). But if there is variation between your replicates (e.g. lab mice, people, etc.), the Poisson assumption will tend to underestimate the variance and any differences you observe (e.g. testing the null hypothesis that two groups have the same mean expression) will be overstated.

If DE is not your interest, then you'll have to explain what your end goals are (at least to me).

HTH,
Mark (co-author of edgeR)
Question: Hi there. I am running delly to detect CNV.

There are 6 BAM I need to run. Is there any difference if I run it like

delly -t DUP -o sample1.out -g ref.fa sample1.bam
delly -t DUP -o sample2.out -g ref.fa sample2.bam
...
delly -t DUP -o sample6.out -g ref.fa sample6.bam

# and then merge the output vcfs

Versus

delly -t DUP -o allSamples.out -g ref.fa sample1.bam sample2.bam sample3.bam sample4.bam sample5.bam sample6.bam

Thanks

Answer: Second sentence of the methods section in the Delly paper:

> All input libraries are analyzed jointly to achieve optimal sensitivity.

specificity > sensitivity : call separately

sensitivity > specificity : joint call
Question: For many this will be very basic, but for me it's driving me nuts!

I have a tab-delimited header:
```
sampleA	sampleB	sampleC
```

Depending on the number of samples I have, I want to make a new header as such:

For 3 samples:
```
sampleA	sampleB	sampleC	output001	output002	output003
```

For 4 samples:
```
sampleA	sampleB	sampleC	sampleD output001	output002	output003 output004
```
and so on...

I can duplicate the initial header to get the correct number of columns:
`paste header header`

However, I want new column names.. and that's where I'm stuck....

Answer: Use `awk` and `NF`. You can use `printf` to pad 0s.

Example:

echo -e "sampleA\tsampleB\tsampleC" | awk -F"\t" -vOFS="\t" '{printf($0"\t"); for(i=1;i
Question: I invested some time into learning R and would prefer not to learn Python or Ruby just for convenience of parsing XML.

I am not a pro R programmer so I need a higher level library then XML package of R.
getNodeSet and XPath is very hard to use.

Answer: raygozak answered the first part of your question: Bioconductor is the nearest equivalent for R of the Bio* libraries. Unfortunately, it does not include a package to handle EUtils queries.

I appreciate your difficulties when using R to parse XML; compared with other languages it is not very intuitive and the documentation is poor (lots of it but rather technical). However, I think you may just have to bite the bullet and get to grips with the R XML library. There are a few good tutorials and blog posts around on the Web (see comments under your question).

There are also some good examples in Chapter 8 of the book Bioinformatics and computational biology solutions using R and Bioconductor, some of which you can read at Google Books.
Question: Hello everyone,

I have the following fastq files:

- Sample_1_R1.fastq.gz
- Sample_1_R2.fastq.gz
- Sample_2_R1.fastq.gz
- Sample_2_R2.fastq.gz
- Sample_3_R1.fastq.gz
- Sample_3_R2.fastq.gz
- Sample_4_R1.fastq.gz
- Sample_4_R2.fastq.gz
- Sample_5_R1.fastq.gz
- Sample_5_R2.fastq.gz

I'm wondering if anyone know a perl or bash script that allows me to run BWA-mem simultaneously on those samples. My output has to be:

- Sample_1.sam
- Sample_2.sam
- Sample_3.sam
- Sample_4.sam
- Sample_5.sam

Answer: Here is a basic function that is then called in parallel for each sample using `GNU parallel`. Using `-j` in `parallel` allows to decide how many jobs run in parallel. Please see the manual of `parallel` for details. The output will be a BAM file for each sample (there is no need/advantage keeping SAM files as they are not compressed and only take up space).

function BWA {

INDEX=$1
BASENAME=$2
bwa mem "${INDEX}" "${BASENAME}"_R1.fastq.gz "${BASENAME}"_R2.fastq.gz | samtools view -o "${BASENAME}"_aligned.bam

}; export -f BWA

ls *_R1.fastq.gz | awk -F "_R1.fastq.gz" '{print $1}' | parallel "BWA /path/to/bwa/index {}"
Question: Hi!

I have ped/map file containing two populations (~110 samples, ~500k SNPs). Is there a software which can phase the data? Alleles should be coded as 1 and 0, where 0 reflects ancestral allele.

I am trying to detect selection and at the moment preparing an input file for selscan programm. So if you can suggest good alternative programs it will help as well.

Thanks.

Answer: Was the data in VCF format before plink? There are lots of ways to phase in VCF. BEAGLE is my preferred phasing algorithm.

Also if your are interested, I have been writing selection software that works directly from VCF.

https://github.com/jewmanchue/vcflib/wiki
Question: How can I find the relation between pairs of pathways in pathway databases like KEGG?

Is there any R-package to do this?

Answer: I believe that the needed data is available through KGML files. For example the [hsa05130 pathway][1] has few arrows pointing on other pathways, for example [hsa04210][2]. If you'd get the 5130 as a [KGML file][3] through their REST API, you'll find that these arrows are encoded as:

```

```

i.e. as *external links*. Should not be difficult to parse an XML using R. As a quick solution:

```
require(XML)
hsa05130_kgml <- xmlParse("http://rest.kegg.jp/get/hsa05130/kgml")
out_nodes = getNodeSet(data, "//entry[@type='map']", fun=xmlToList)
links_summary <- data.frame(do.call(rbind, out_nodes))
head(links_summary)
```

```
1 Regulation of actin cytoskeleton, #000000, #FFFFFF, roundrectangle, 805, 460, 119, 34
2 Flagellar assembly, #000000, #FFFFFF, roundrectangle, 126, 253, 113, 25
3 Adherens junction, #000000, #FFFFFF, roundrectangle, 761, 975, 120, 25
4 Tight junction, #000000, #FFFFFF, roundrectangle, 631, 975, 98, 25
5 Toll-like receptor signaling pathway, #000000, #FFFFFF, roundrectangle, 584, 264, 132, 34
6 Bacterial secretion system, #000000, #FFFFFF, roundrectangle, 271, 361, 185, 34
.attrs
1 1, path:hsa04810, map, http://www.kegg.jp/dbget-bin/www_bget?hsa04810
2 4, path:map02040, map, http://www.kegg.jp/dbget-bin/www_bget?map02040
3 5, path:hsa04520, map, http://www.kegg.jp/dbget-bin/www_bget?hsa04520
4 6, path:hsa04530, map, http://www.kegg.jp/dbget-bin/www_bget?hsa04530
5 7, path:hsa04620, map, http://www.kegg.jp/dbget-bin/www_bget?hsa04620
6 8, path:map03070, map, http://www.kegg.jp/dbget-bin/www_bget?map03070
```

[1]: http://www.genome.jp/kegg-bin/show_pathway?hsa05130
[2]: http://www.genome.jp/kegg-bin/show_pathway?hsa04210
[3]: http://rest.kegg.jp/get/hsa05130/kgml
Question: Hi,

Imagine that I have 20 bam files, which I want to change their header then merge them and then sort them! So I just wondering if these two scenarios make different output :

1. reheader (each of them) -> merge all -> sort the output
2. merge all -> reheader the output -> sort

All bam files have the same read group, and I want to add sample name to each RG field

@RG ID:1120

Thanks

Answer: It shouldn't make any difference, but what I'd suggest doing is using `samtools merge` and passing the header you want to apply to the merged file through the `-h` argument. That will simplify your workflow. `samtools merge` does expect your input to be already sorted though. If your input is indeed sorted, then that will allow you to accomplish your goal with a single command.

For example, put the header you want in a file called `header.sam`:

samtools merge -h header.sam merged.bam 1.bam 2.bam 3.bam ... 20.bam