cannin/biostars_qa huggingface dataset:
Question: Hey guys!

I did a multiple sequence alignment with sequences from proteins with disordered regions (IDR). My goal is to see if there is conservation between the 20 proteins for the disordered regions. I already have all the information concerning those IDR since I did the disorder prediction myself. My problem is that I would like to highlight the specific amino acids that are disordered in Jalview using the data that I already have. I do not want to select the regions manually so I was wondering if there is a way to import annotations using my data.

I looked in the manual and I found the basics of it all but I could not find anything that would let me do this automatically, not with the cursor.

I hope this explains my problem clearly and I'll try to find a solution as soon as I can. If someone has an easy answer that would be awesome! :)

Have a good day!

Answer: It all depends what format your disorder predictions are in. If you can massage them into a [Jalview annotations file][1] like so:

```
JALVIEW_ANNOTATION
# Created: Tue May 26 10:53:09 BST 2015
SEQUENCE_REF UniRef50_Q9C8G6
LINE_GRAPH IUPredWS (Long) Protein Disorder with IUPredWS - raw scores
Above 0.5 indicates disorder 0.5707,0.5707,[72144e]|0.5807,0.5807,[72144e]|0.4619,0.4619,[72144e]|0.4685,0.4685,[72144e]|0.4801,0.4801
GRAPHLINE IUPredWS (Long) 0.5 Above 0.5 indicates disorder ff0000
```

Then you can just load it onto your alignment, and then use the Colour->By annotation dialog with the 'per sequence' tick box selected, and 'Above Threshold' on the dropdown menu ([see the documentation for a figure][2]).

[1]: http://www.jalview.org/help/html/features/annotationsFormat.html
[2]: http://www.jalview.org/help/html/colourSchemes/annotationColouring.html
Question: I am working on results of CNV studies obtained form an online database. I have the start and end co-ordinates of CNVs of several individuals. I want to visualize and compare all of them chromosome wise and identify the CNV regions common in all the individuals. How do I go about this? Is there an online software which would allow this? PS: I am not a trained bioinformatician.

Answer: You can also attach this data to the [Ensembl Browser][1] as explained in this [tutorial page][2]. You will have your data in a format like 'chr start end individual' (i.e. [BED format][3]) and you can compare it to our [gene annotation][4], our [variation][5] data (variants annotated according to [consequence terms][6]), and [regulation data][7] from ENCODE, Blueprint and Roadmap projects (for annotation of enhancer promoter, CTCF, etc). We also include the CNVs from 1000 Genomes and structural variants from DECIPHER (see the region where *[PRR29][8]* is anntoated). Once you attach your data, you also have the option to share that visualisation with your colleagues and collaborators. You just send them the URL we generate and they will see your data against our annotation.

[1]: http://www.ensembl.org
[2]: http://www.ensembl.org/info/website/tutorials/userdata.html
[3]: http://www.ensembl.org/info/website/upload/bed.html
[4]: http://www.ensembl.org/info/genome/genebuild/index.html
[5]: http://www.ensembl.org/info/genome/variation/index.html
[6]: http://www.ensembl.org/info/genome/variation/predicted_data.html#consequences
[7]: http://www.ensembl.org/info/genome/funcgen/index.html
[8]: http://uswest.ensembl.org/Homo_sapiens/Share/8f0b56753b831f8614c104d1a780f16611452918
Question: Hi,

Does anyone come across this error in `picard SortSam` command? ([site](https://broadinstitute.github.io/picard/index.html#GettingHelp))

The command run:

java -Djava.io.tmpdir=$TMP_SPACE -jar $PICARD \
SortSam I=$IN_BAM O=$OUT_BAM SORT_ORDER=coordinate \
VALIDATION_STRINGENCY=LENIENT

The command tries to sort a bam file by coordinate. It runs for 2 hours and then it shows the error below:

> [Thu Jan 07 00:14:51 UTC 2021] picard.sam.SortSam done. Elapsed time: 150.85 minutes.
Runtime.totalMemory()=24461180928
To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp
Exception in thread "main" java.lang.IllegalArgumentException: Value (4135) to large to be written as ubyte.
at htsjdk.samtools.util.BinaryCodec.writeUByte(BinaryCodec.java:331)
at htsjdk.samtools.BAMRecordCodec.encode(BAMRecordCodec.java:155)
at htsjdk.samtools.BAMRecordCodec.encode(BAMRecordCodec.java:40)
at htsjdk.samtools.util.SortingCollection.spillToDisk(SortingCollection.java:254)
at htsjdk.samtools.util.SortingCollection.add(SortingCollection.java:182)
at htsjdk.samtools.SAMFileWriterImpl.addAlignment(SAMFileWriterImpl.java:187)
at picard.sam.SortSam.doWork(SortSam.java:161)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:304)
at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:103)
at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:113)

I search a bit, but I'm not a java programmer. As far as I can tell, it does not seem to be a problem of space neither RAM memory (I run this command in one server with big memory). The file in question is quite big (a few hundred GB).

The error seems related with a problem to write into the disk, but I'm not familiar with what is `ubyte` neither how to interpret and possibly solve the problem (if possible).

Thank you in advance for help or advice,

António

Answer: > htsjdk.samtools.BAMRecordCodec.encode(BAMRecordCodec.java:155) at

the error is here: https://github.com/samtools/htsjdk/blob/master/src/main/java/htsjdk/samtools/BAMRecordCodec.java#L155

it's thrown when the read NAME is too long (the SAM spec says length(name) < 256 (coded as uint8_t) )

so it's a problem with your bam. Check the length of your read names.
Question: Dear Biostars,
I want to find a way to search for reads that have reverse complements inside a library from a RIP-seq experiment.
I have found this related question: [Question: Finding complementarity between mRNA and rRNA using BLAST][1]
but nucleotide blast seems to be a bit different now.
Current blastn searches for reverse complement alignments too?
Is there a package that I could do it in R?

thank you in advance

[1]: https://www.biostars.org/p/197164/

Answer: One option is `bedtools intersect` or `intersectBed` which can be run on two bam files.

intersectBed -a shortReads.bam -b longReads.bam -S -bed

The -S flag will force the **opposite** strandedness (as opposed to the `-s` flag that forces same strandedness). You'll want to play with the `-f` flag that requires the fraction of A that must be overlapped (the default is 1bp, you may want to raise `-f 1.0` to require the entire short read overlaps the long read).
Question: In Gene ontology and KEGG analysis, Why there are two cutoff values; p-value and (-) log10(FDR) value?

What is the role of these two cut-offs?

Could bioinformatician give an explanation to a wet lab person?

Answer: I'm going to link to [another of my previous answers](https://www.biostars.org/p/352631/#352666) that explains the logic behind enrichment analyses. The reason that both p-value and false discovery rates (FDR) are used is to account for statistical issues arising from [multiple testing](https://en.wikipedia.org/wiki/False_discovery_rate). In short, p-values are computed from a single score, which can be misleading due to the high number of tests actually performed.

Multiple test correction methods, FDR among them, try to correct for this issue to give more accurate values that can be compared to whatever threshold (or **alpha**) that you choose - typically 0.05 or 0.01. The chosen threshold is arbitrary, but those values are generally regarded as stringent enough to be relatively believable. A score with a p-value or q-value or FDR of 0.05 calculated would occur in 5% of cases by random chance in the data. Using lower alpha values should reduce the false positive rate with the drawback of increasing the false negative rate. [This short paper](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2907892/) has a simple example that is fairly easy to follow along with additional context.
Question: My lab has sequenced an exome using the "Agilent SureSelect All Exon Kits" ( http://www.genomics.agilent.com/CollectionSubpage.aspx?PageType=Product&SubPageType=ProductDetail&PageID=3041 )

Does anyone knows where I can find a BED describing the captured regions ?

P.

Answer: Ok https://earray.chem.agilent.com/suredesign/index.htm was the right place to find the BED. I was able to log-in today

Thank you all.
Question: I came across an [article][1] benchmarking different tissue-specific gene expression identifying methods.

I am wondering why is it not a common practice to use differential expression analysis methods, like DESeq2, to identify tissue specific genes? Practically, it can be really a time-consuming task, if you have multiple tissues - identifying DE genes using one-vs-one for each tissue and then find the common DEs for the tissue of interest. Apart from practicality, is there any reason not to do it?

[1]: https://doi.org/10.1093/bib/bbw008

Answer: Philosophically, tissue *specific* suggests a gene is _only_ expressed in one tissue. Just because a gene is DE between tissue type A and B, with positive logFoldChange (i.e. higher in tissue A), that doesn't mean it is not expressed in B, it can easily be expressed in both, while being higher in tissue A.

Practically, most DGE tools assume that:

1. That changes between the two conditions are limited, or at least the mean change is 0.
2. That the compositional changes more or less reflect the expression level changes. The normalisation methods in DESeq2 and edgeR are designed to cope with a certain amount of compositional change between samples (i.e. changes in the total amount of RNA in a cell between samples/changes in expression levels of highly expressed genes changing the reads available for other genes), but there is probably a limit to how big the changes can be.
3. That the "gene" is the same sequence in both samples - i.e. there haven't really been any big splicing changes (although using tximport and weighting data be effect gene length can offset this somewhat).

All three of these assumptions are less likely to be true when comparing tissues than when comparing two samples from the same tissue.
Question: I am trying to align the following sequence using a python code. I am able to obtain the best score results but not the alignments with lower scores. Is there an efficient way to do this?

# Import pairwise2 module
from Bio import pairwise2

# Import format_alignment method
from Bio.pairwise2 import format_alignment

# Define two sequences to be aligned
X = "ACAGCCAGCGGGGGAAATAGAGTATTTTTTTCGAGGAGGATAGGTTACTCCTCTGGTCAAAGTTGGATGCGTAATATCACACGATGCGTTGACCGCGGAGAGTGACATTATGTGGATCGGAATATCACTATTCTGCTCGAACTTCCATAG"
Y = "ACAGCCAGCGGGGGAAATAGAGTATTTTTTTCGAGGAGGATAGGTTACTCCTCTGGTCAAAGTTGGATGCGTAATATCACACGATGCGTTGACCGCGGAGAGTGACATTATGTGAATATCACTATTCTGCTCGAACTTCCATAGAGATCG"

# Get a list of the global alignments between the two sequences ACGGGT and ACG satisfying the given scoring
# A match score is the score of identical chars, else mismatch score.
# Same open and extend gap penalties for both sequences.
alignments = pairwise2.align.globalms(X, Y, 1, -1, -1, -0.1, force_generic=True, score_only = True)

# Use format_alignment method to format the alignments in the list
for a in alignments:
print(format_alignment(*a))

Here is the result I get with the above code.

ACAGCCAGCGGGGGAAATAGAGTATTTTTTTCGAGGAGGATAGGTTACTCCTCTGGTCAAAGTTGGATGCGTAATATCACACGATGCGTTGACCGCGGAGAGTGACATTATGTGGATCGGAATATCACTATTCTGCTCGAACTTCCATA------G
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| |||||||||||||||||||||||||||||| |
ACAGCCAGCGGGGGAAATAGAGTATTTTTTTCGAGGAGGATAGGTTACTCCTCTGGTCAAAGTTGGATGCGTAATATCACACGATGCGTTGACCGCGGAGAGTGACATTATGT------GAATATCACTATTCTGCTCGAACTTCCATAGAGATCG
Score=141

What I want is alignments with lower scores such as 130, 100, etc. I know there could be many such examples but I need the alignments with the lower scores. I tried looking at the source code for Biopython but it seems to only give the best score.

Answer: At the end of [**this paper**][1] there is a short program for global alignment by dynamic programming. If you know how to modify C code, it may help in your experiments.

> What I want is different scores for the same match, mismatch and gap penalties.

There is only one global alignment for the same match, mismatch and gap penalties. If you want to get multiple alignments from identical parameters, you will have to go with local alignments - and even in that case there may be only one. These algorithms are meant to go through all possible alignments, but the shortest backward path requirement means that they will output only one alignment. Look through the code I referenced above, and you may be able to modify it so it outputs a fixed number of alignments according to decreasing scores.

[1]: https://www.nature.com/articles/nbt0704-909
Question: Hi guys I'm in the following situation: My dataset: 400 patients: 150 of them had surgery between 1994-1998 and 250 had surgery between 1998-2005. Of the 150, 50 had leukemia, 30 breast cancer and 70 had lung cancer. Of the 250, 65 had leukaemia, 120 breast cancer and 65 lung cancer. The question is ti check if there's an enrichment of a tumour type according to the range of years of surgery (1994-1998 and 1998-2005). I was thinking to perform a Fisher test but, for some reasons, my boss would like a logistic regression. I have no idea how to perform the logistic regression with such data. Can anyone help me please?

Answer: If you're only interested in whether there's a difference in proportion of the two types of cancers between the two groups, then a chi-squared test or Fisher's exact test are appropriate and simpler. A logistic regression is possible but it would evaluate the cancer type as a function of the year, i.e. it would answer whether the year predicts the type of cancer, which is not really what you seem to be asking from the data. If you don't know why you should do one approach over the other, I suggest you ask your supervisor for explanation, maybe there are more factors that need to be taken into account.
Anyway, you can do logistic regression with the glm() function in R, something like

model <- glm(cbind(breast,lung) ~ year, data = your.data.set, family = binomial)
Question: Hi all,

Any script available for doing sequence length disrtibution of fastq file?

Thanks,
Deepthi

Answer: Here is a solution using `awk`:

awk 'NR%4 == 2 {lengths[length($0)]++} END {for (l in lengths) {print l, lengths[l]}}' file.fastq

It reads like this: every second line in every group of 4 lines (the sequence line), measure the length of the sequence and increment the array cell corresponding to that length. When all lines have been read, loop over the array to print its content. In `awk`, arrays and array cells are initialized when they are called for the first time, no need to initialize them before.
Question: Hi
I am very new to bioinformatics analysis. I am trying to map a whole genome in galaxy, but I don't know exactly how to find my reference genome (chlamydomonas reinhardtii). I do know It is available on NCBI but don't know which file should be used as a reference genome. I mean: is the reference genome only one single file, or is it multiple files for each chromosome?
I'd be grateful if you could help me.

Answer: You can find the representative genome page for your organism at NCBI ([LINK][1]). Actual genome sequence is [in this file][2] which is in fasta format. This is the file you should use as a reference.

If you need specific help with Galaxy then consider posting questions to their help forum: https://help.galaxyproject.org/

[1]: https://www.ncbi.nlm.nih.gov/genome/147?genome_assembly_id=22881
[2]: https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/002/595/GCF_000002595.1_v3.0/GCF_000002595.1_v3.0_genomic.fna.gz
Question: Given a directory containing multiple FASTQ files, I would like to retrieve the run information from each file (e.g. flowcell, lane number, index e.t.c) and export the data to a CSV file. Before I write a script to do this myself, is anyone aware of any software which does this already?

Answer: I wrote **illuminadir** for my lab http://lindenb.github.io/jvarkit/IlluminaDirectory.html . output is a json or xml file;

```
$ find . -type f -name "*.fastq.gz" | java -jar dist/illuminadir.jar -j
{
"samples": [
{
"files": [
{
"forward": {
"file-size": 82852,
"md5filename": "3854678b2381caa92758a980eb0d180d",
"path": "./src/test/resources/SAMPLE1_GATGAATC_L002_R1_001.fastq.gz",
"side": 1
},
"id": "p1",
"index": "GATGAATC",
"lane": "2",
"md5pair": "6fb0b6d5dec57436724e8efcb736aed5",
"reverse": {
"file-size": 82625,
"md5filename": "5450dea531fb2d807402d6bcc59cf15d",
"path": "./src/test/resources/SAMPLE1_GATGAATC_L002_R2_001.fastq.gz",
"side": 2
},
"split": "1"
}
],
"sample": "SAMPLE1"
}
],
"undetermined": []
}
```
Question: Hello,

I am new to ATAC-analysis, and wish to ask for some advice in identifying significant differences between two ATAC-seq samples.

Firstly, I have two conditions - control and knockdown.
I have mapped the PE reads to a reference genome, and then subsequently taken reads lower than 100bp (i.e. open chromatin regions).
I then use MACS2 to turn my two un-normalised BAM files to peaks using this command:

macs2 callpeak -t macs2 callpeak -t bamfile --outdir /path/to/ -f BAMPE --keep-dup all --pvalue 1e-2 --call-summits --bdg

I then use feature counts to quantify mapped sequencing reads my 'peak' genomic features.

awk 'OFS="\t" {print $1"."$2+1"."$3, $1, $2+1, $3, "."}' peaks.narrowPeak > contol.peaks.saf

featureCounts -a peaks.saf -F SAF -T -p -o control.peaks_countMatrix.txt control.lessthan100bp.bam

So I result in two matrix outputs from feature counts: control.peaks_countMatrix.txt; and knockdown.peaks_countMatrix.txt

I am not used to using DESeq2 for differential analysis (have been using Kallisto/Sleuth for RNA-seq). Am I right in thinking that the two matrix files can be inputted to DESeq2, to peform differential analysis against? Does anyone have a pipeline that they could share in getting this done?

Answer: I would call peaks with the command:

macs2 callpeak -t bamfile --nomodel --keep-dup=all

The `--nomodel` turns off the ChIPseq-specific shifting model. Do not use the p-value argument, as correction for multiple testing prevents you from excessive numbers of false-positives. This excessive number of peaks will contain lots of regions with small counts (because of the relaxed cutoff) which lack the statistical power to call differential accessability. Better use the default q-value cutoff of 0.05 or even 0.01. As you already have the sub-100bp reads, I think the `--call-summits` option is neither necessary, nor beneficial. Most peaks will already be around 100-200bp.

There are multiple options how to call the peaks. Either call them one the merged BAM files of one condition, or call on every replicate and then take the intersect for each condition. There is no standard for this. I would probably call peaks with the above command and a q-value of 0.01 on each single file (I hope you have replicates, if not there is no point in using DESeq2), and then take the intersect between the files of one condition. You'll end up with one peak set for each condition, which is then to be merged (`bedtools merge`) and used as input for featureCounts to count the reads. DESeq2 takes **one** count matrix which contains the counts for each single replicate over one peak set (the merge between the conditions). From there on, you can pretty much follow the [DESeq2 manual][1]. It contains all the necessary code and is outstandingly extensive and informative.

[1]: https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
Question: I want to know if there is any free open source tool available for Gasteiger-Hückel charge calculation? I know it is available on Sylbyl but is there any other(free) way to calculate these charges for small molecules?

I was actually trying this paper, if any one is interested,

How to improve docking accuracy of AutoDock4.2: a case study using different electrostatic potentials.

Hou X1, Du J, Zhang J, Du L, Fang H, Li M.X et al.

Answer: Did you have a look at [Acpype][1], there is an option (`-c` gas) that should calculate these charges on a ligand. You only need to download Antechamber (part of the Ambertools package) and acpype.

[1]: http://svn.code.sf.net/p/ccpn/code/branches/stable/ccpn/python/acpype/
Question: Does anyone know of any software package that can convert .impute2 data to .mldose (i.e. imputed data from IMPUTE2 to imputed data from MACH)? I have tried impute2mach in GenABEL, but it has repeatedly failed with a known error...

Direct conversion preferred, but I'm happy with indirect as long as it works!

Answer: IMPUTE output is posterior probabilities, 3 columns for 1 SNP, e.g.: `0.1 0.1 0.8` corresponding to `AA AB BB`. Meaning it is more likely to be `BB`.

MACH outputs dosage, 1 column for 1 SNP, e.g.: `1.7`, corresponding `0=AA, 1=AB, 2=BB`, so `1.7` is more likely to be `BB`.

Depending on the size of the data, you can write a quick script in R that would convert posterior probabilities to dosage:

0.1*0 + 0.1*1 + 0.8*2 = 1.7
Question: Hi all

I have multiple csv files stored in a list (list.data) having same column names but different data frame names.
I would like to subset a column (CHR) by a specific value (1) and apply it to all the data frames in the list and to all the columns.

I wrote this code :

CHR1<-lapply(list.data, function(x) { x["CHR"] = 1; x }

However this seem to subset correctly that specific column (CHR) but the other columns did not change, subset was not applied to the rest of the columns.

I wrote also this code but it returns an empty list

CHR1<-lapply(list.data, subset, "CHR"=1)

Any help highly appreciated.
Thanks

Answer: Try:

CHR1 <- lapply(list.data, function(x) { x[ x$CHR == 1, ] })

- `=` is an assignment operator, e.g.: `myVar = 1; print(myVar)` same as `<-`.

- `==` is a logical operator, used for comparing, e.g.: `1 == 2` this will give us `FALSE`.
Question: Hello everyone,

I am trying to use Varscan to find somatic variants but I want to apply filters on the minimum variant frequency and minimum number of variant reads. I have tumor-normal paired samples.

In this case, should I use:
a) mpileup2snp and mpileup2indel separately and indicate the filters
or
b) somatic and then somaticFilter on each one of the files produced (snps and indels)?

Is there a difference in how these two approaches work?
Also, I want the output in vcf format and according to the manual, only mpileup2snp/mpileup2indel offes that option?

Thanks!

Answer: It is the second workflow that you want. I have [a script][1] in my Github that uses VarScan2 for somatic calling (not saying it is a nice one or fulfills any good-practice standard when it comes to style or whatever), but you might get some inspiration from it. It also uses `GNU parallel` to parallelize the process over all chromosomes.

It starts by calling the raw variants with `VarScan2 somatic`, then separates them into somatic and germline with `processSomatic`, selects high-confidence variants with on `VarScan`'s Fisher's Exact Test, and finally runs the heuristic `fpfilter` to remove junk calls.

[1]: https://github.com/ATpoint/Variant_Calling/blob/master/VarScan2/VarScan2_Somatic.sh
Question: Hi,

Could you explain briefly what's the goal of the "dereplication" process in a metagenomics analysis ?

Answer: I think I understood the purpose !
Derepliaction is not a filtering process ! This step is used to compute the abundance of reads.
For instance, I have the following fasta file :

> seq1
ACGT
> seq2
ACGT
> seq3
TCGA
> seq4
GGAC

After dereplication I will have :

> seq1_2
ACGT
> seq3_1
TCGA
> seq4_1
GGAC

Note that abondance value are written in the name of sequence