cannin/biostars_qa huggingface dataset:
Question: I have two bacterial srains - one is wild type and the other has one gene deleted. In my experiment, I calculate some phenotype for both of these strains. I repeat this experiment 5 times. Now I want to see if there is a difference in this phenotype between the two strains, so I perform a t-test (assuming the data is normally distributed for now). My question is - should I perform a paired t-test or unpaired? I am a bit confused about it since I feel I have some rationale for doing it either ways -

Why I think its **paired** : in every replicate, I compare the same two strains and I repeat this 5 times

Why I think its **unpaired** : Its not like I am subjecting the same strain to two different conditions (which often happens in paired tests). Rather, I have two different strains being subjected to the same condition.

Any help would be appreciated!

Answer: Assuming your use of the t-test is valid, you should use the unpaired version. Your question is getting at the difference between biological and technical replicates, which is different than paired vs. unpaired measurements.

A paired measurement would be taken across conditions or time (the groups being compared), so for example measuring 5 individuals at two time points and testing the distributions of time1 vs. time2. This is different than what you're doing, which is doing technical replication within each group but comparing those to independent observations in another group (5 replicates in group1 versus 5 replicates in group2).

I think you reasoned through this already in your statement "I think its unpaired."

You would, however, expect there to be less variation under your technically replicated experimental design than if you had used 5 biological replicates. You might mention this in your limitations.
Question: Hello I'm trying to run bowtie for a several files using a loop so I don't have to repeat the command over and over. The loop looks something like this:

mm10="GenomeFolder/Bowtie2Index/genome"

for (( i = 76; i <= 79; i++ ))
do
bowtie2 -x $mm10 -1 SRR51363$i_1.fastq.gz -2 SRR51363$i_2.fastq.gz \
done

The problem is that I can't use the $i variable in between SRR51363 and _1.fastq.gz, It causes an error. Any adivce on the right syntax? I tried to google this but it's hard to describe this problem I'm having. Thank you!

Answer: Put curly braces around variable names, and quotes around words:

```
mm10="GenomeFolder/Bowtie2Index/genome"

for (( i = 76; i <= 79; i++ ))
do
bowtie2 -x "${mm10}" -1 "SRR51363${i}_1.fastq.gz" -2 "SRR51363${i}_2.fastq.gz"
done
```

Curly braces in this context prevent variable names from being expanded incorrectly. Quotes allow files and paths to contain spaces and special characters. While the latter may not be an issue here for this specific example, it is a good habit to get into.
Question: Hi,

How can I remove duplicate rows from the output of my toptableF? I have tried the following without success:

if (unique(fit2$genes)){
Sigenes.F <- topTableF(fit2, number=100, genelist=fit2$genes, adjust.method="bonferroni",
sort.by="F")
}

Thanks

Answer: You'll need to pay attention to what you are doing. If you are new to R, you will need to do some reading on R data structures, as each behaves differently with regard to indexing. Just guessing, but you probably want:

SigGenes.F[!duplicated(SigGenes.F$genes),]
Question: I have some doubts regarding my understanding of blast algorithm that i'd like to dispel. It's well known that blast tries to find the same words in query and database and then uses it as starting point for the alignment extension. Let's imagine we have word AGTcCAT in query sequence and one word AGTgCAT in database. So there is a single mismatch between the words. My question, may the blast use this pair as starting point for the alignment extension? I think not, since blast requires only perfect matches of two words. Am i right? Thanks!

Answer: AFIAK, the seed region can't have mismatches. And the seed is created from the query string, which is matched against the database to find potential hits, which are then extended on both sides to get the final alignment. In your case, if the word size of seed region is >4, then you won't be able to find any hits. If seed <=3, there are matches (AGT = AGT; CAT=CAT etc), which can be extended according to mismatch penalty.

Note: For nucleotide, the defulat word-size is 11 and for AA, it is 3.
Question: I've aligned 454 RNA-Seq reads using bwa. Does somebody know if a software that identify variants that suffered splicing exists ?

Answer: TopHat is probably not the "best" in terms of sensitivity or specificity. Have a look here for a more complete list. I use GSNAP quite regularly for this task.
Question: When running `samtools idxstats` on an indexed BAM file `foo.bam` in order to build a list of chromosomes, I get a result that includes a line I'm not sure how to interpret:

```
$ samtools idxstats foo.bam | cut -f1
chr1
chr10
chr11
...
chrY
*
```

What is the asterisk supposed to refer to, when I iterate through this list and extract data from that chromosome?

Answer: in the code of `bam_index.c`

printf("*\t0\t0\t%llu\n", (long long)idx->n_no_coor);

where `n_no_coor` = "unmapped reads without coordinate"
Question: Hi, I have a .vcf containing 1 million SNPs for a subset of a population (80 individuals) and what I would like to do is select SNPs that are informative for building pedigree and then redesign a probe set to target those SNPs for the rest of the population.

I am not quite sure how to determine which SNPs out of the 1 mil are actually kinship informative. I do know parent-offspring and full/half sib relationships for this subset of individuals but will need to be able to build a pedigree out to 4th degree relatives for the entire population. Does anyone have any suggestions on a good approach for selecting these SNPs? Would it be something similar to how I would select ancestry informative (which I have already mined out of this SNP dataset)?

Thanks!

Answer: I'd suggest LD pruning/clumping to retain the unlinked SNPs (PLINK?), followed by identifying IBD (Identity-By-Descent) using PLINK. Pruning essentially finds representative SNPs in LD blocks to capture genotype information for the entire blocks and reduces the amount of variant data to analyze.

https://privefl.github.io/bigsnpr/articles/pruning-vs-clumping.html

https://www.cog-genomics.org/plink/1.9/ld#indep

https://www.cog-genomics.org/plink2/ibd

This would result in relatedness (kinship) measures. See this:

https://www.biostars.org/p/58663/

But sometimes the relatedness measure could be a bit different from pedigree/parentage relationships you're after, as relatedness measures are also affected by population histories like genetic drift in small populations and/or bottlenecks.
Question: Hello!

I'm super noob on DESeq interactions so sorry in advance if this question is duplicated. I just need to ask directly to see if I understand correctly.

I have two experiments performed in young and old mice. I want to test in both for knock down vs control of a certain gene. I want to make a differential expression analysis between conditions (ko vs ctrl) and between experiments (young vs old).

My coldata looks like

sample condition experiment
mouse1 ctrl young
mouse2 ctrl young
mouse3 ctrl young
mouse4 ko young
mouse5 ko young
mouse6 ko young
mouse7 ctrl adult
mouse8 ctrl adult
mouse9 ctrl adult
mouse10 ko adult
mouse11 ko adult
mouse12 ko adult

But I am not really interested on what are all the differences between young and old mice, but **only those that have an impact on the ko vs ctrl test**.

So my design for the DESeq object is ~ condition + experiment + condition:experiment

> resultsNames(dds)
[1] "Intercept" "condition_ko_vs_ctrl" "experiment_young_vs_adult" "conditionko.experimentyoung"

I am not sure, but the intuitive thing to do would be to take experiment_young_vs_adult. Is this correct?

Also, could someone explain why do I not have "conditionko.experimentadult", "conditionctrl.experimentadult", "conditionko.experimentyoung" in resultsNames(dds)?

Thank you in advance!

Answer: I found this website helpful in working out the right way to deal with interactions

https://rpubs.com/ge600/deseq2

I think you likely want "conditionko.experimentyoung", or what's likely more useful is to relevel age so that young is the control, and then the contrast will be "conditionko.experimentadult"
Question: Hi,

I've got some 450K data and I've been looking at the various packages around....There's a lot of normalisation methods.

I realise normalisation as a method is quite subjective to the data, and it's probably best to "shop around", and see which method best fits my data. However, are there any methods that in principle should be better than others?

Thanks,

Answer: Hi Andrew,

you are corrected, depending on your dataset characteristics, some normalization methods are more appropriate than others. If the goal of your analysis is to compare different tissues, cell types, or cancers versus normals (or any other dataset with global methylation differences), you might want to take a look at our recent paper comparing normalization methods for such datasets: http://www.genomebiology.com/2014/15/12/503. We showed that Functional normalization performs really well at conserving global methylation differences (`preprocessFunnorm()` function in the minfi package)

Otherwise, if you are doing a EWAS-type of study, you may want to use `preprocessNoob()` + `preprocessQuantile()`, both implemented in the minfi packages as well. The former performs a background correction based on the negative control probes (implementation of the noob method described in http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3627582/) that is always appropriate to use (and btw is called internally by `preprocessFunnorm` as well). The `preprocessQuantile()` function is an adaptation of quantile normalization to the 450k array, very appropriate when no global differences are expected between your samples (for instance whole-blood samples).

Finally, if you have a large dataset and think you could have batch effects, there is a recent paper that applies a batch effect correction method (RUVm) in top of Funnorm: http://nar.oxfordjournals.org/content/early/2015/05/18/nar.gkv526.full

Hope this helps,

Jean-Philippe
Question: I have been trying to run OMA standalone using a few genomes downloaded from the OMA browser and then adding a couple of my custom genomes. All the steps of OMA standalone including All-vs-All computations, pairwise orthologs etc have successfully finished. However, I have consistently faced the same error with the final step of generating and storing HOGs for these genomes.

Here is the error message that I get despite trying to run OMA on different machines:

********************************************************************************

Storing results
Traceback (most recent call last):
File "/home/orthologs/OMA.2.4.1/bin/..//.venv/bin/warthogs.py", line 128, in
run_gethogs()
File "/home/orthologs/OMA.2.4.1/bin/..//.venv/bin/warthogs.py", line 95, in run_gethogs
Settings.check_consistency_argument()
File "/home/orthologs/OMA.2.4.1/.venv/lib64/python3.6/site-packages/gethogs/settings.py", line 53, in check_consistency_argument
cls.inputfile_handler = file_manager.inputfile_handler_factory()
File "/home/orthologs/OMA.2.4.1/.venv/lib64/python3.6/site-packages/gethogs/file_manager.py", line 97, in inputfile_handler_factory
return OmaStandaloneFiles(os.path.normpath(join(settings.Settings.pairwise_folder, '..')))
File "/home/orthologs/OMA.2.4.1/.venv/lib64/python3.6/site-packages/gethogs/file_manager.py", line 44, in __init__
for row in csv_reader:
_csv.Error: field larger than field limit (131072)
********************************************************************************
An error occured in bottom-up HOG computations:

--- WARTHOGs:
- Start at 13:36 on 2020-12-08
- Orthology relations folder:Output/PairwiseOrthologs (standalone format)
- Method use to merge HOGS: pair
- Output file: Output/HierarchicalGroups.orthoxml

Error, (in GetHOGsBottomUp) failed to compute bottom-up hogs

Can somebody please help me with this error?

Also is there a way to only run the HOG computation stage of the program without having to generate all the output files again?

Answer: Hi Gaurav, (hi Alberto)

I had a look at the problem and it turns out that there are some proteins with a huge number of crossreferences in the fasta header. The csv module has a limit on the length of each field, which is exceeded by these proteins.
The fix is quite simple, so you can actually do it yourself if you want. I also plan to release a new version of oma standalone within the next days that contains that fix (will be 2.4.2).

what you need to do is to include in `hog_bottom_up/gethogs/file_manager.py` right after the import lines:

csv.field_size_limit(1<<20)

After this change, the analysis run through in a few minutes without further problems.

Thanks for pointing out the problem and good look with our analysis

Adrian
Question: Hi,

The below is my Dmrcate run and I do get the DMRS. but how to get the CPGids that are part of these DMRS

corfit <- duplicateCorrelation(myMs, design=designMatrix, block = biolrep)
dmrCate.<-cpg.annotate("array", myMs, what="M", arraytype = "EPIC", analysis.type="differential", design=designMatrix, coef=2, block=biolrep, correlation=corfit$consensus.correlation,fdr=0.01)
dmrcoutput_all <- dmrcate(dmrCate,lambda=500, C=5, pcutoff="fdr",min.cpgs=1,mc.cores=1)
wgbs.ranges_all.snp7 <- extractRanges(dmrcoutput_all, genome = "hg19")
write.csv (wgbs.ranges_all.snp7,file="DMRCate_DifferentailRange.txt")
write.csv(dmrcoutput_all.snp7$results,file="DMRCate_ALL.txt")

how to get the 17 Cpgids from the result shown below

"","seqnames","start","end","width","strand","no.cpgs","minfdr","Stouffer","maxbetafc","meanbetafc","overlapping.promoters"
"1","chr6",31539539,31540750,1212,"*",17,2.06654851019045e-120,2.16758469969414e-149,-0.162022961391297,-0.106211449278892,"LTA-001, LTA-002, LTA-004, LTA-003"

I tried with gRanges but unable to understand how to perform it

P.S: I have posted the same in Bioconductor Forum. sorry for redundancy. Any help soon is appriciated

Answer: There would be multiple ways(more efficient way ) to what do than the below script which I have used solve my own problem.

But It would be great to know If there is any efficient way to do the same.

**OBJECT** <- is the last object from which beta/mvalue for DMRCate is give as input

gr <- granges(OBJECT)

mydata<-read.csv(file="DMRCate_Output.csv",header=TRUE,stringsAsFactors = FALSE)
for(n in 1:nrow(mydata))
{
chr<-mydata2[n,2]
start=mydata2[n,3]
end=mydata2[n,4]
infos<-names(gr[seqnames(gr) == chr & start(gr)>=start & end(gr) <=end ])
allcpg<-toString(infos)
mydata$cpgids<-allcpg
}

write.table(mydata,file="DMRCate_Model3_DifferentailRange_V3V4_withCpgids",sep="\t");

P.S: Here the output from DMRcate is written to a file
Question: Hi All,

Any solution to plot multiple `heatmap` in one page in R (`heatmap.2` in `gplots`)? `par(mfrow=c(2,2))` don't work. It will plot the heatmap in separate pages.

Thanks.

Answer: I find the solution eventually, Anyway, share with others, hope it is helpful in your routine analysis. Also, I found `gridGraphics` and `grid` is very useful for complicated figure preparation.

library(gridGraphics)
library(grid)

grab_grob <- function(){
grid.echo()
grid.grab()
}

data<-read.table("heatmap.input.matrix.data.txt")

library(gplots)
gl <- lapply(1:ncol(data), function(i){
heatmap.2(data.matrix(data[[i]]))
grab_grob()
})

grid.newpage()
library(gridExtra)
grid.arrange(grobs=gl, ncol=5, clip=TRUE)
Question: Hello Everyone,

I am looking for the best way to get linkage information from unphased whole genome population data. I have a vcf file with multiple individuals from different populations. The data is unphased but I would like to detect regions with an excess of linkage disequilibrium as a measure of positive selection. I have not phased the data because I have a limited number of individuals per population of a non-model species and therefore worry that phasing will be very inaccurate.

What do people think would be the best way to detect regions with high levels of linkage disequilibrium? I was thinking something like VCFtools `--geno-r2` option might be suitable.

Thanks for your help!

Best regards,
Rubal

Answer: The --geno-r2 option in vcftools should be enough for your needs; however, you can not calculate linkage disequilibrium if your data is not phased. If you were studying human individuals, I would suggest you to impute more genotypes by merging it with the 1000 Genomes data, but if you say that you are working with a non-model organism, this is not an option. Is there a close model organism that you can use to impute data?
Question: Hi,

I was wondering if any one could help me. I am currently working on *C. difficile* and have obtained a LuxS mutant that does not always behave as one would expect. I have therefor sequenced it together with the wild type strain to see what is going on. This is about as far as I have got. I used both breseq, bowtie2 and samtools to check my sequences against the published sequence data on ncbi. As however, the *luxS* gene was knocked out through the insertion of an ~600bp insert I have been told that this method of analysis may not be able to determine whether this insert is present or where it is.

If anyone has any ideas on how I would go about determining whether this insert is present within the *luxS* gene. I would be most grateful.

Answer: I would do this:

1. If you know the sequence of the insert, concatenate it to the reference.
2. Map the reads to the modified reference, with an aligner that produces global alignments with high sensitivity (like BBMap), or that produces chimeric alignments (like bwa-mem).
3. Look at the gene in IGV. If there is a long insertion event, you will see some place in it where the reads hit a junction. With global alignments, they will look error-free up to some point then suddenly all the bases will not match the reference. That point is where the insertion occurred, and the part of the read that does not match the reference should match the insert.
4. If your reads are paired, you should see improper pairs such that one read mapped to the gene of interest, and the other read mapped to the insert (which is why it was added to the reference).
Question: Has anyone here worked with CNCI before? I'm just about exhausted trying to figure out what I'm doing wrong.

So, I tried the following:

python CNCI.py candidate_lncs.gtf -g -o test -m ve -p 16 -d ./dbase GRCm38.primary_assembly.genome.fa

and I recieved the following error:

name@server %57> CNCI.py -f candidate_lncs.gtf -g -o test -m ve -p 16 -d ~/dataBase/GRCm38.2bit
/home/pgian1/src/CNCI
sh: 1: twoBitToFa: not found
Traceback (most recent call last):
File "CNCI.py", line 533, in
GtfInFiles = open(fastaFiles)
IOError: [Errno 2] No such file or directory: 'candidate_lncs.gtf.fa'

It looks to me that CNCI is trying to read my GTF file as a fasta file. I'm not sure why though. I have tried different versions of python, changing my path, etc and nothing is working. Any help at all is appreciated!! Or if someone could find a protocol or something that demonstrates CNCIs usage. I have found nothing so far

Answer: CNCI is trying to open your files as a Fasta file, becauase the first thing the program does if you provide GTF files is convert them to fasta. Almost certainly this conversion is failing and CNCI is therefore unable to open the FASTA file it thinks it should have created.

Looking at the CNCI source code, the conversion is done in two steps. The first uses a perl script from CNCI's own directory to convert the GTF file to a bed file. The second uses twoBitToFa to convert said bed file to a FASTA file. The failure could be either one of these two steps.

Thus I would try to the following commands to see which fail:

perl PATH_TO_CNCI//gtf2Bed.pl candidate_lncs.gtf > candidate_lncs.gtf.bed
twoBitToFa -bed=candidate_lncs.gtf.bed ~/dataBase/GRCm38.2bit candidate_lncs.gtf.fa
Question: Hi Guys,

I want to concatenate all these values in the columns enclosed with "" and separated by comma in the order as shown below. How can I get this done in R?

`df1`

```
chr start end strand
chr4 443333 232444 +
chr5 4455332 4433323 -
chr5 4443333 4433355 +
```

I want this result

Result (exactly as shown below):

"chr4","443333","232444","+","chr5","4455332","4433323","-","chr5","4443333","4433355","+"

Answer: You have to use the `t` (transpose) function.

```
> paste(dQuote(t(df1)), collapse=', ')
""chr4"," 443333"," 232444","+","chr5","4455332","4433323","-","chr5","4443333","4433355","+""
```
Question: Hi,

I finished mapping and variant calling on these for solid data. Now I'm looking to identify large rearrangements and transposable elements. Can anyone tell me any good tools for these? One tool I have come across now is nFuse. I wanted to know what the general pipeline for this is.

And I would like to view such arrangements as Circos plots.

Thanks for the help.

Answer: There are a number of non-reference TE insertion mapping tools available now including:

- [RelocaTE](https://github.com/stajichlab/RelocaTE)
- [ngs_te_mapper](https://github.com/bergmanlab/ngs_te_mapper)
- [PoPoolationTE](https://code.google.com/p/popoolationte/)
- [TE-locate](http://sourceforge.net/projects/te-locate/)
- [RetroSeq](https://github.com/tk2/RetroSeq)
- [TEMP](https://github.com/JialiUMassWengLab/TEMP)

These systems differ in how they map NGS reads, how they discover and represent non-reference TEs in a genome (see https://www.biostars.org/p/49563/), and in their final output file formats. We have tried to characterise the major features of each of these systems (besides TEMP) in this poster. There has been no systematic evaluation of all of these methods, so it is difficult to say which is the "best".

In the absence of a thorough evaluation study, our view is that it is necessary to run multiple TE mapping tools and evaluate the impact that different tools have on your biological conclusions. To this end, we are currently developing a meta-pipeline called McClintock that installs and runs all of the tools listed above (besides TEMP) and outputs results in a standard BED format.