cannin/biostars_qa huggingface dataset:
Question: Hello everyone! This is my first time posting here so I hope I am doing this the right way.

I'm a student in bioinformatics in Québec, Canada and I just started a summer project with my teacher. I have to predict disordered proteins from a database containing a lot of sequences. I don't know if any of you is familiar with this process but I'll ask anyway.

I'll use DISOPRED3 for the prediction. The problem is that the results I get are not the same on my personal computer as those I get from the server itself. Before using my scripts to predict with all the files I have, I would like to have the same results as the server. I use UniRef90 (like they do) and I installed the last version of the program. What could be different and cause the small differences? I already sent an email asking what they are using but they didn't respond yet.

Thank you, any help will be appreciated!

Answer: If you read the FAQ at the psipred download page (http://bioinfadmin.cs.ucl.ac.uk/downloads/psipred/), you get:

"2. The next most common question we get regards getting different predictions from our PSIPRED web server from those you get on your own system. There are many reasons why this might occur, but the most obvious one is that we are using a
​different sequence data bank to the one you are using. Because modern secondary structure prediction methods are based on

analysing multiple sequence alignments, if your data bank includes some extra sequences or misses out some sequences compared
to our local data bank then you might get a slightly different prediction. Hopefully the differences will be small, but
they can be quite significant if the alignment only includes a few sequences. So, if our server can find say just 5
homologous sequences and your system finds 20 homologous sequences to align, then the predictions may be very different
indeed. That's just the reality of analysing evolutionary information. If the alignments are different, then the predicted
secondary structure will probably be different.

It's also possible that we are running a slightly older version of PSIPRED. Oddly enough we don't update our servers immediately
after releasing a new version of PSIPRED as we need to do internal testing first. So, check that you are running the exact
same version of PSIPRED as our server is currently running. Even so, it's likely that the real problem is going to be down
to differences in the data banks and alignments that are produced."

Also, I am no expert with this software but some structure modelling software are not always deterministic and have a constant seed parameter for debugging purpose. Others (like Rosetta Antibody) are wonderful but will never run or even compile correctly on your laptop.
Question: Dear all,
I want to do GATK variant calling. I have several raw-paired-end-reads fastq files for one pig individual. And I have finished mapping using bwa mem. GATK best practices (for calling germline SNPs + indels) tells that after mapping, I need to merge bam files. And it recommends `picard MergeBamAlignment`. what's the difference between `picard MergeBamAlignment` and `samtools merge` please? I found the input files they need are different. The former one need unmapped bam file as well as the reference genome file. And in my understanding, each `MergeBamAlignment` command is only for one fastq (pair) file, not for multiple fastq files. Is there any other difference please? Particularly, are their aims, and their output results any different please? I feel that `MergeBamAlignment` is for combining the information both from mapped file and unmapped file. samtools merge is for combining diffferent mapped file into one big mapped file. Or Can I ignore GATK's advice and simply use `samtools merge` instead? Thank you.

Yingzi

Answer: Hi Yingzi,so you are here!

I haven't used `MergeBamAlignment(Picard)` but I think the best way to get knowledge about a software is to read the original documentation.

As for your case, firstly, we can find the asscoiated description in the documentation of [MergeBamAlignment(Picard)](https://software.broadinstitute.org/gatk/documentation/tooldocs/4.0.0.0/picard_sam_MergeBamAlignment.php) and the documentation reads:
>A command-line tool for merging BAM/SAM alignment info from a third-party aligner with the data in an unmapped BAM file, producing a third BAM file that has alignment data (from the aligner) and all the remaining data from the unmapped BAM. Quick note: this is not a tool for taking multiple sam files and creating a bigger file by merging them. For that use-case, see {@link MergeSamFiles}.

Secondly, we can also find the documentation of [samtools](http://www.htslib.org/doc/samtools.html) and we can find the function of merge:
>Merge multiple sorted alignment files, producing a single sorted output file that contains all the input records and maintains the existing sort order.

As you can see,the two tools have different functions, so you need to use both of them. What's more, you can also merge multiple bam files using [MergeSamFiles(Picard)](https://software.broadinstitute.org/gatk/documentation/tooldocs/4.0.0.0/picard_sam_MergeSamFiles.php) instead of `samtools`. You can also find that the desctiption of `MergeSamFiles(Picard)`:
>Merges multiple SAM and/or BAM files into a single file. This tool is used for combining SAM and/or BAM files from different runs or read groups into a single file, similarl to the "merge" function of Samtools (http://www.htslib.org/doc/samtools.html).

So I can confidently to say that what you have understood is right.

Yours :)

Duo
Question: Hi everyone,

I have exome data of 25 unrelated patients and 40 unrelated control samples. I'm looking for rare variants associated with the disease.

I've checked some obvious things: one variant present in cases not in controls. One gene enriched for rare variants in cases compared to controls.

However, now I want to do some statistical testing to find association of a variant/gene with the disease. Since this is just a small-scale study I don't think I can work with methods used in GWAS. Can someone point me to some papers/methods/ideas that could be of interest for my specific situation?

Thanks in advance

Answer: As Sam pointed out, SKAT is a good place to start but unless there is a **very** strong signal from one gene/variant it's unlikely that you will get a statistically significant association from 25 cases. [This paper][1] from Lander & co has some calculations to estimate minimum sample size needed in these types of studies to get good results, and even under the best conditions there generally needs to be hundreds of cases. But it's worth a shot and even if nothing comes up with a low enough p-value the top hits still might be of interest.

If you haven't already, try answering questions such as: "what rare variants show up in *at least two* cases, but no controls", "what genes contain novel variants in *at least two* cases, but no controls"

Have you done any filtering by types of variants (eg, removing all synonymous variants from your analysis, or only looking at stop loss/gain variants)?

[1]: http://www.pnas.org/content/early/2014/01/16/1322563111
Question: Hi all,

When I'm using `FeaturePlot()` for one gene, my graph has the color gradient guide but when I want to plot multi genes, the color gradient guide is missed.

FeaturePlot(seurat.HSC, features ='Ecm1', min.cutoff = 'q10')

![enter image description here][1]

vs

FeaturePlot(seurat.HSC, features = c('Ecm1','Gucy1b1','Acta2'), min.cutoff = 'q10',
split.by = 'sample' )

![enter image description here][2]

How to solve it? Thanks for any help.

[1]: /media/images/9bbc9437-c384-4bc8-a874-3ee0b01f
[2]: /media/images/a0ae54a0-658a-4dc7-bfb7-da0079e4

Answer: Hello,

I think the problem is related to *split.by*, try to use this:

FeaturePlot(seurat.HSC, features = c('Ecm1','Gucy1b1','Acta2'), min.cutoff = 'q10', split.by = "sample") & theme(legend.position = "right")

Best
Question: I am attempting to use Salmon (version 0.12.0) to quantify transcript counts from RNAseq data (unstranded paired-end reads).

The existing pipeline in my lab is to trim/qc fastq files, align them with STAR (2.5.2a), merge sample BAMs together (if they were spread across lanes), sort them with samtools, and then feed them to Salmon in alignment-based mode. STAR was run with the `--quantMode TranscriptomeSAM` option and the `genomeDir` pointed to a genome generated using STAR's `genomeGenerate` function with the `--sjdbGTFfile` option pointing to a GTF file.

Here's the Salmon run line:

salmon quant -t /hg38_salmon_transcriptome.fa -l IU -p 16 -a some.transcriptome.sorted.bam -o ./

This produces a quant.sf file that seems to make sense, although it also produces a massive (15+ GB) error file that seems really upset about suspicious pairs (see here for someone else with the same issue: https://www.biostars.org/p/164823/ )

WARNING: Detected suspicious pair ---
The names are different:
read1 : K00274:68:HGYCHBBXX:5:1119:13311:37220
read2 : K00274:68:HGYCHBBXX:3:2103:3325:33598

I was a bit sick of this behavior and decided to give Salmon the fastq files directly in alignment-free mode.
This pipeline was to trim/qc my fastqs, combine the files from different lanes and then gives the reads to salmon. The program runs fine and produces no errors.

Here's that Salmon run line:

salmon quant -i /hg38_salmon_transcriptome_index -l IU -1 samp_R1.fq.gz -2 samp_R2.fq.gz -o ./

The big issue is that the quant.sf files look really different between these two pipelines. Like the transcript ENST00000361739 (this is the MT-CO1 gene) has a TPM of 40735 in alignment-free mode, but a TPM of 105 in alignment-mode. Further, the range of TPMs (and counts) varies widely between the two modes. In alignment-free mode, I'm getting TPMs in the thousands, but in alignment mode the highest is 300 and most values are under 100.

So, my questions are:

1. Does anyone know what's causes the "suspicious pair" error when I'm running in alignment-based mode?
2. **Why am I getting such huge differences in the outputs of these two strategies?**

Answer: The command line you show for running using the STAR alignments:

```
salmon quant -t /hg38_salmon_transcriptome.fa -l IU -p 16 -a some.transcriptome.sorted.bam -o ./
```

suggests that the bam file provided to salmon was (coordinate?) sorted; is this the case? In alignment mode, salmon, like e.g. RSEM, expects all of the alignments for a read to be adjacent in the BAM file, and for the mates of a pair to follow each other. STAR will output the alignments like this by default, as long as you don't ask it to sort the output. Could you verify if this is the case and, if so, see how things look if you don't coordinate sort the output bam?
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: all <- paste(apply(df1, 1, function(row) paste(dQuote(row), collapse=",")), collapse=',"+",')
Question: Hello, Everyone,

The raw data is coming from https://satijalab.org/seurat/v3.1/pbmc3k_tutorial.html.

pbmc3k_filtered_gene_bc_matrices.tar.gz.

It is the single-cell file: barcodes.tsv, genes.tsv, matrix.mtx.

I wanted to transfer matrix.mtx file to csv file. I tried many python mat2csv, none of them worked.

Thanks in advance for any help!

Best，

Yue

Answer: R

library(Matrix)
matrix_dir = "/home/li/"
barcode.path <- paste0(matrix_dir, "barcodes.tsv")
features.path <- paste0(matrix_dir, "features.tsv")
matrix.path <- paste0(matrix_dir, "matrix.mtx")
mat <- readMM(file = matrix.path)
feature.names = read.delim(features.path,
header = FALSE,
stringsAsFactors = FALSE)
barcode.names = read.delim(barcode.path,
header = FALSE,
stringsAsFactors = FALSE)
colnames(mat) = barcode.names$V1
rownames(mat) = feature.names$V1

Python

import csv
import gzip
import os
import scipy.io
matrix_dir = "/home/li/"
mat = scipy.io.mmread(os.path.join(matrix_dir, "matrix.mtx"))
features_path = os.path.join(matrix_dir, "features.tsv")
feature_ids = [row[0] for row in csv.reader(open(features_path), delimiter="\t")]
gene_names = [row[1] for row in csv.reader(open(features_path), delimiter="\t")]
feature_types = [row[2] for row in csv.reader(open(features_path), delimiter="\t")]
barcodes_path = os.path.join(matrix_dir, "barcodes.tsv")
barcodes = [row[0] for row in csv.reader(open(barcodes_path), delimiter="\t")]
Question: Hi,

SAM file specification has fields/tags such as alignment score (AS:i: in bowtie2) and edit-distance (NM:i:) for each end of the pair.

I am working with amplicons paired-end data (thus overlap is not uncommon) and would like to know if there is a paired-end "combined" score for the whole mate/pair (taking into account the overlap), e.g. for alignment score, edit distance, etc; this will be useful for later filtering etc.

I'm currently using bowtie2, and have just noticed "The alignment score for a paired-end alignment equals the sum of the alignment scores of the individual mates." written in the manual. Is that indeed the case especially when there is an overlap? How are the scores calculated for each end (0.5 for each ?)

I'm currently using my own implementation as a post-processing which is inefficient.

many thanks
Shim

Answer: Overlaps are not in any way taken into account when calculating the `AS` tag in bowtie2. There is no score reported that's a direct sum of the mate's `AS` score, though one can sort of use the MAPQ as an ersatz value for this, since MAPQ is dependent on the summed `AS` from the mates (as well as their summed `XS` values). But really if you want to compute a joint `AS` and such while accounting for overlaps then post-processing is your only option.
Question: Hello, I have a list of human genes and I'd like to retrieve the physical coordinates (GRCh37/hg19 assembly) of their 3'UTRs. Are you aware of any software that can do that? Thanks !

Answer: For RefSeq annotation, you can use the `add_utrs_to_gff` python script to first add 5' and 3' UTR features and then use unix `grep` to extract the genes of your interest. The latest RefSeq annotation of the GRCh37 assembly is here: https://ftp.ncbi.nlm.nih.gov/genomes/all/annotation_releases/9606/105.20190906/

## download annotation in GFF3 format
$ curl -O https://ftp.ncbi.nlm.nih.gov/genomes/all/annotation_releases/9606/105.20190906/GCF_000001405.25_GRCh37.p13/GCF_000001405.25_GRCh37.p13_genomic.gff.gz

## download the add_utrs_to_gff3 python script
$ curl -O https://ftp.ncbi.nlm.nih.gov/genomes/TOOLS/add_utrs_to_gff/add_utrs_to_gff.py

## add utr features to the gff3 file
$ python3 add_utrs_to_gff.py GCF_000001405.25_GRCh37.p13_genomic.gff.gz > GRCh37_with_utrs.gff3

## extract 5' UTR for GeneID:5768
$ grep 'five_prime_UTR' GRCh37_with_utrs.gff3 | grep -w 'GeneID:5768'
NC_000001.10 BestRefSeq five_prime_UTR 180123968 180124042 . + . ID=utr00100412821;Parent=rna-NM_001004128.2;transcript_id=NM_001004128.2;Dbxref=GeneID:5768,Genbank:NM_001004128.2,HGNC:HGNC:9756,MIM:603120
NC_000001.10 BestRefSeq five_prime_UTR 180124004 180124042 . + . ID=utr00282651;Parent=rna-NM_002826.5;transcript_id=NM_002826.5;Dbxref=GeneID:5768,Genbank:NM_002826.5,HGNC:HGNC:9756,MIM:603120
Question: Hello,

I am trying to analyze the public dataset https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE126030
I've downloaded the fastq files onto my cluster, and would like to proceed with cellranger count.

I am in a test folder and the only file is: SRR8526547_1.fastq and refdata-cellranger-GRCh38-1.2.0

cellranger count --id=cellranger \
--transcriptome=/home/jl2/scratch60/refdata-cellranger-GRCh38-1.2.0/ \
--fastqs=.\
--sample=SRR8526547_1.fastq \

I keep getting the error of

Invalid path/prefix combination: /gpfs/ycga/scratch60/k/jl2/test, ['SRR8526547_1.fastq']
No input FASTQs were found for the requested parameters.

Can't seem to figure out what's wrong. Does it need fastq.gz instead of fastq?

Answer: Hello,

I have been troubleshooting
> error: No input FASTQs were found for the requested parameters.

for several hours now. In my case the file names, the file path and the command were all fine. Finally what solved the issue for me was to move the fastq.gz files into a seperate folder that only contained fastq.gz files. The original folder had some other files in it (md5, fastqc output, etc.). Not sure why this was a problem for the pipeline, but make sure to give this a try if you run into similar trouble.

Best,
Max
Question: Hello,

I'd like to study how NCBI's non-redundant protein database (nr) has developed over the years. However, I'm yet to find a way to download anything but the latest release from the NCBI ftp. Are those old versions lost for good from the public domain?

Answer: As far as I am aware NCBI do not provide archived versions of the 'nr' database, although they might be available upon request.

However since most of the sequences in 'nr' come from the protein translations in GenBank and UniProt provide archived releases for UniProtKB (which includes translations from EMBL-Bank), the UniProt releases would probably cover what you need. See ftp://ftp.uniprot.org/pub/databases/uniprot/previous_releases/.

Alternativly the UniProt's UniParc database is equivalent to the NCBI's 'nr' database, and provides additonal date information which would allow you to create subsets based on the database at a particular date. For the XML version of the UniParc database, which contains the additional information, see ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/uniparc/

Please note: the NCBI 'nr' database and the UniParc database are sets of non-identical sequences (i.e. the database contain one sequence for each unique sequence, with meta-data providing details of all the source entries containing the sequence). Non-redundant sequence databases such as UniRef or those generated with CD-HIT are different, and merge subsequences such as those from sequencing fragments into either the longest or a representitive sequence. To generate your own 'nr' like database(s) use the 'nrdb' program (http://blast.advbiocomp.com/pub/nrdb/) on your collection of sequences.
Question: Why do we remove duplicates from BAM files while using Samtools?
When we have paired end data we can remove duplicates as a fragment OR as a pair. How do each of these methods differ?

Answer: I would personally recommend using Picard for marking or removing duplicates. If you have a paired data, then both reads for a pair will be used to select duplicates. In this case, if there is another pair that has both of its reads aligning at the same exact location as this pair, then one of these would be marked as duplicates. For fragment reads, location of only one read will be used to mark the duplicates.
Question: Hello !

I don't know what the difference between " ENSG00000002586.1" and "ENSG00000002586.19_PAR_Y" in Ensembl gene/transcripts ID means.

How do you deal with these ID with "_PAR_Y" in RNA-seq mapping process ?

Thanks,

Answer: It means the gene has multiple copies in the [pseudoautosomal regions][1]. The '.1' part is the version number.
How you deal with them most likely depends on whether you care about these regions.

[1]: https://www.ensembl.org/info/genome/genebuild/human_PARS.html
Question: Hi,

I have two lists of genes

mycounts <- read.csv("geneID and length.csv", header = T, sep = "\t", stringsAsFactors = FALSE)

```
> colnames(mycounts)
[1] "genesID.geneslength"

> head(mycounts[1:4,])
[1] "R0010W,1272" "R0020C,1122" "R0030W,546" "R0040C,891"

> dim(mycounts)
[1] 7130 1
```

mycounts1 <- read.table("read.txt", header = T, sep = "\t", stringsAsFactors = FALSE)

```
> dim(mycounts1)
[1] 5961 1
> colnames(mycounts1)
[1] "Freq"
```

How I can have only genes in my read file in my genes file? I mean genes file has 7130 that I only need 5961 of them

May you help me please? Thank you

Answer: Something is weird here, since you have the names files joined with its length and separated by a comma. Am I right?. I believe that since you have only 1 column as I can see from the dim(). It is likely that you need to do the `read.csv` in a different way to be able to separate both values, like using a different sep value. I need to know the format of the original file to suggest you

In that case none of these suggestions will work because the genes name are common, but then the length value has to be the same to do a merge
Question: I have `VCF` file (containing diallelic variants) and reference genome in `Fasta` of some non-model plant. I'd like to extract `SNPs` flanking sequences. I've found that `bedtools` and `samtools faidx` could be to some extent useful, but apparently don't solve the issue. I need to get output (`fasta` or tabular file) with following SNPs representation :
TCTCTGCCAATCACTAGAGGCCGCTTTCGCTTTTA[A/G]TTTGTGTGTGGTCAGAGTTCTTCCGGACTTT

Answer: Here a little python solution:

import pysam

# open vcf file
vcf = pysam.VariantFile("input.vcf")
# open fasta file
genome = pysam.FastaFile("genome.fa")
# define by how many bases the variant should be flanked
flank = 50

# iterate over each variant
for record in vcf:
# extract sequence
#
# The start position is calculated by subtract the number of bases
# given by 'flank' from the variant position. The position in the vcf file
# is 1-based. pysam's fetch() expected 0-base coordinate. That's why we
# need to subtract on more base.
#
# The end position is calculated by adding the number of bases
# given by 'flank' to the variant position. We also need to add the length
# of the REF value and subtract again 1 due to the 0-based/1-based thing.
#
# Now we have the complete sequence like this:
# [number of bases given by flank]+REF+[number of bases given by flank]
seq = genome.fetch(record.chrom, record.pos-1-flank, record.pos-1+len(record.ref)+flank)

# print out tab seperated columns:
# CRHOM, POS, REF, ALT, flanking sequencing with variant given in the format '[REF/ALT]'
print(
record.chrom,
record.pos,
record.ref,
record.alts[0],
'{}[{}/{}]{}'.format(seq[:flank], record.ref, record.alts[0], seq[flank+len(record.ref):]),
sep="\t"
)

fin swimmer

**EDIT**:

I rewrote the code. The version before had an off-by-one error.
Question: I'd like to get my hands on a microarray dataset with as many conditions as possible. One thought that I had was to simply concatenate many individual microarray datasets together, however it seems that they often dont have the same probes, and sometimes are in different units (fold change, intensity, etc).

Is there a way I can download all the datasets of GEO, but only include datasets that have the same probes, and are normalized in the same way?

Answer: Did you have a look at the Expression Atlas? It integrates many gene expression experiments to make them comparable.
Question: Hi,

I'm a medical student with genetics, and I have come to understand some knowledge in bioinformatics would help me getting a research position in the field. I've been browsing the internet for some course, but all the ones I've found so far did not suit my needs. Does anyone know of an online bioinformatics course that:

1. is affordable
2. is open to an undergraduate medical student
3. provides a certificate
4. teaches you the most important things for genetic, genomic and molecular biology research

Answer: The [Data Analysis for Life Sciences series from EdX][1] meets all of your criteria. I've been too busy to enroll, but I did get my hands on some of the notes and lectures and they're top-notch.

[1]: https://www.edx.org/course/data-analysis-life-sciences-1-statistics-harvardx-ph525-1x