cannin/biostars_qa huggingface dataset:
Question: I'm involved in establishing a new bioinformatics lab. It will serve primarily for teaching purporses, but I think it is possible to design small grid as well.
The idea is to have 15-20 desktop computers, 4-8 cores each, and build grid using half of the cores. The budget is ~30k Euros. Probably, some network file system (NFS) will be needed as well.
Can you share yours experience in that matter? Or maybe you have some suggestions of ready solutions? Or is it better to buy cheaper desktops and invest saved money into 2-3 very strong workstations?

What I experienced so far in Bioinformatics labs, can be divided into 2 solutions:

all data stored in network file system (NFS) and computer clients loading everything from NFS. Independent cluster uses the same NFS. It has advantage of easy back-ups, and unified system. But the users have limited privileges and have to ask for every piece of software. Unfortunately, in case of troubles with NFS, no one can work at all:/
the client computers store data locally, but there is NFS, mainly for cluster and back-up purposes

I'm awaiting yours comments.

Answer: Populating 15+ separate desktops with a number of bioinfo packages does not seem right to me. While it is doable, you would need something like cfengine or Chef to automate everything.

You have many options when it comes to exporting particular directories from NFS server, from everything (remote booting of clients) through /biosoft and /home dirs, or just /biosoft. So some work can be done during NFS failure assuming you got the basic tools in your /home dir.

Given a choice I would go with dumb and cheap clients (but with enough RAM and swap & /tmp) and few servers with a lot of RAM / cores. Some stuff (i.e. genome/transcriptome assembly) hardly works, if at all, on low RAM machines.

If you are concerned about NFS failures then you may go with RAID and mirrored NFS servers.
Question: Hi, I have a reference genome (fasta) and list of genotypes (vcf) for one individual. I want to make a haploid genome file by randomly selecting one of the individual's alleles at each variable site. I would use GATK's alternate reference maker, but I do not want to bias my choice to either the reference or non-reference allele. Anybody know how I can do this?

Thank you!

Answer: ```
#!/usr/bin/perl

use strict;
use warnings;
use autodie;

$ARGV[2] or die "use newGenome.pl \n";

my $vcf = shift @ARGV;
my $fas = shift @ARGV;
my $out = shift @ARGV;

# load sequences
my %seq;
$/ = "\n>";
open (my $fh, "<", $fas);
while (<$fh>) {
s/>//g;
my ($id, @seq) = split (/\n/, $_);
$seq{$id} = join "", @seq;
}
close $fh;

# parse VCF
$/ = "\n";
open (my $vh, "<", $vcf);
while (<$vh>) {
next if (m/^#/);
chomp;
my ($chr, $pos, $id, $ref, $alt, $qual, $fil, $info) = split (/\t/, $_);
next unless ($alt =~ m/^[ACGT]$/); # only valid SNPs
my $afrq = 0.5; # in case of missing values use a 50% chance
if ($info =~ m/AF=(0\.\d+)/) {
$afrq = $1;
}

my $prob = rand();
if ($prob < $afrq) {
substr ($seq{$chr}, $pos - 1, 1) = $alt);
warn "modified $chr:$pos $ref -> $alt\n";
}
}
close $vh;

# Write modified sequences
open (my $oh, ">", $out);
while (my ($id, $seq) = each %seq) {
print $oh ">$id\n";
while ($seq) {
print $oh substr ($seq, 0, 80), "\n";
substr ($seq, 0, 80) = '';
}
}
close $oh;
```

\* Cold night, procastinating work, perl lover
Question: Hey

can u help me with getting only one representative fasta sequence per family?
Is there way to simply do that?

cheers

X

Answer: It's not trivial. You could use the sequences from the trRosetta Pfam model set, which are representative of the family ([download link][1]).

We have a method for getting representative sequences in [our paper][2] if you are comfortable with using hmmsearch:

> A representative target sequence was found for each family using hmmsearch to search the UniRef90 database with the Pfam HMM and taking the closest subsequence match by E-value.

[1]: http://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/structure_models
[2]: https://www.nature.com/articles/s41467-019-11994-0
Question: I have R1.fastq and R2.fastq from RNA-Seq experiment. I'm currently trying to perform differential expression using Salmon and DESeq2.

Since the RNA-Seq was done using stranded library kit, it was under my assumption that R1.fastq contains forward reads and R2.fastq contains reverse reads.

So, I used the following command to run Salmon:

salmon quant -i -l ISF -1 -2 -o --gcBias

This resulted in poor mapping rate (~3%) and confused me at first since I saw no issues with the sequencing data (e.g. no rRNA contamination, etc...).

Instead, when I tried with below command, I was able to reach ~90% mapping rate:

salmon quant -i -l ISR -1 -2 -o --gcBias

From the documentation, F/R of option ***-l*** tells Salmon which strand a read is coming from given a .fastq file:

> F: read1 comes from the forward strand and read2 comes from the reverse strand
>
> R: read1 comes from the reverse strand and read2 comes from the forward strand

So does this mean that in my R1.fastq, there are reverse strand reads?

Thanks in advance.

Answer: > Since the RNA-Seq was done using stranded library kit, it was under my
> assumption that R1.fastq contains forward reads and R2.fastq contains
> reverse reads.

That assumption was wrong, as you can see with your own data. In Truseq standed kits, read 1 is in the reverse orientation.
Question: Hi,

I have a sample with a strange name. It has been giving me some issues with some software that I use. I am trying to change the sample name in the VCF file using bcftools reheader

The command I am using is

bcftools reheader -s "99686 C99686\n" input.vcf -o output.reheader.vcf

but when I try looking at the samples name in the file using bcftools query, the name is unchanged.

bcftools query -l /gpfs_backup/meurs_data/TestBCF/output.reheader.vcf

I assuming I am doing something wrong with the reheader command but I am not sure what I am doing wrong.

Answer: the option -s of bcftools reheader is a **FILE** not a string. https://samtools.github.io/bcftools/bcftools.html

-s, --samples FILE
new sample names, one name per line, in the same order as they appear in the VCF file. Alternatively, only samples which need to be renamed can be listed as "old_name new_name\n" pairs separated by whitespaces, each on a separate line. If a sample name contains spaces, the spaces can be escaped using the backslash character, for example "Not\ a\ good\ sample\ name".
Question: I am trying to retrieve [pLI score][1] for some human genes but I am not able to find a DB of this score. Does anybody knows how to retrieve such data for a specific gene?

[1]: https://genome.ucsc.edu/cgi-bin/hgTrackUi?db=hg19&g=gnomadPLI

Answer: The pLI score (now deprecated by the ExAC o/e score). You can pull these tracks (pLI and oe) from UCSC:

https://genome.ucsc.edu/cgi-bin/hgTrackUi?db=hg19&g=gnomadPLI

and downloaded here:

http://hgdownload.soe.ucsc.edu/gbdb/hg19/gnomAD/pLI/
Question: Hi everyone,

I have a wildtype and 2 mutant samples. I have already run lrt (likelihood ratio test) and obtained > 40,000 transcripts but did not get any fold change in the output. After googling, I found that wald test (wt) would be a suitable option. Is it possible to run wald and lrt parrallely or in a comparative manner? What things I should keep in mind while filtering the results from both? How should I proceed further to find out the differentially expressed genes in mutant samples vs. wild?

Thank you!

Answer: Hello there,

If you are doing a simple comparison between wild-type and mutant samples, then you can technically do the wald test. However, in our preprint [1], all of the analyses are done with the likelihood ratio test. This is because we found that the FDR is much more believable under the likelihood ratio test then the Wald test. My general impression is that the Wald test gives too many false positives, although, the rankings are fairly comparable between the two tests.

All that being said, I would suggest simply sticking with the likelihood ratio test.

You do not get any fold change in the output because the test is inherently different. However, you can manually compute the log fold change between the samples because you have a fairly simple design. Let me know if you have any questions about doing this. The function `kallisto_table` should give you a data frame that is easily manipulated with standard tools such as dplyr.

[1] http://biorxiv.org/content/early/2016/06/10/058164.abstract
Question: Hi,

I have a a large dataset which has information in columns of the table of a protein name and a gene name.
I would like to put it online on a public server in some way that is searchable by the scientific community by protein name and by gene name, and possibly links to other established databases such as UniProt.

I saw that there is the possibility of using mySQL, but before diving into learning mySQL (which I am not sure how long it could take), I am curious about what the community of BiosStars suggest, in case there are faster alternatives.

I code in R normally and I saw the R Shiny, could this be a good alternative to mySQL to build an online seacrhable database? What are the pros/cons from people's experiences?
https://cran.r-project.org/web/packages/shiny/index.html

Thank you for your thoughts,
James

Answer: using [DT][1], creating a searchable table on the web,

[shiny][2], building for a webserver,

and [dplyr, interacting with mysql][3], will be a quick way.

[1]: https://rstudio.github.io/DT/
[2]: https://shiny.rstudio.com/
[3]: https://beanumber.github.io/mysql-r-webinar/dplyr.html
Question: I have downloaded the assembly fasta files for all chromosomes for an organism from NCBI. Now I want to map my unique data (WGS) to the genome. How to combine the chromosome-wise fasta files to give a fasta file for the whole genome? This needs to be done to run BWA, right?

Answer: `cat *.fa > genome.fa`

You'll need to index that with bwa or whatever aligner you plan to use then.
Question: Hi,

I am learning goseq and I know a bit of R. I just can read the codes and understand them but not very creative to type what I want to do. So I need to prepare data from deseq2res. One vector containing all genes (column 1) and the other vector containing all DE genes (`padj < 0.05`). Can anyone do the favor? Any lines of code is appreciated!

Answer: d <- read.csv("deseq2res.csv", header=T, row.names=1)
all_genes <- row.names(d)
DE_genes <- all_genes[which(d$padj<0.05)]
Question: I have gone through the deseq2 paper where it tells that its not a simple logfold change of numerator by denominator ,its shrunken LFC which it calculates the maths behind it kind of im trying to understand in simple manner or layman term as i have to give some explanation why rna seq and RTqpcr are not similar...

I have 4 control and 4 test sample respectively all the counts are normalized used deseq2 default

gene baseMean log2FoldChange lfcSE stat pvalue padj	nH4	nH3	nH2	nH1	C4	C3	C2	C1
ENSG00000112773.13	1156.5575015292	2.52671958320643	0.330914786026763	7.63555963619613	2.25E-14	5.31E-12	326.344363070055	316.430980291062	369.993520743799	241.523773415627	1626.54796885324	3322.42939303791	1321.51219870859	1727.67781411328

so how come the log2foldchange or logfoldchange is low.?even though there is a big difference between my sample labelled as H vs C which are stem cells and C's are progenitor cells ?

Another data point ,i have used lfcMLE for logfold change which says "“unshrunken” log2 fold
changes (for a simple comparison, the ratio of the mean normalized counts in the two groups"

> gene baseMean log2FoldChange lfcMLE lfcSE stat pvalue padj SO_7660_BR_02 SO_7660_BR_03 SO_7660_BR_04 SO_7660_BR_06
>
> 797.414 2.972936 3.976479 0.1360428 21.85294 7.288395e-106 1.204553e-101 113.6162 75.41853 1320.512 1680.109

log2foldchange = 2.972936
lfcMLE = 3.976479

there is difference log2foldchange and lfcMLE somehow its still less as compared to qPCR ...

Answer: The fold-change is barely shrunken for these samples, the naive fold-change is ~6.4 and the shrunken fold-change is ~5.8. So I don't know why you think the log2FC is particularly low. The thing with the naive fold-change is that it will be innately unreliable due to the low number of samples and high variability, which is why it's shrunken toward 0 with a normally distributed prior.

The general idea behind shrinking p-values is that fold-changes will generally be overestimated given the actual variance in the samples and the typical miniscule number of samples. So lowly expressed genes are strongly pulled toward 0 while more highly expressed genes are less strongly pulled toward 0.

Regarding a discordance with your qPCR data, without at least knowing what the fold-change was it's impossible for us to judge how actually different the two values are. In general there's good (though not perfect) concordance between the two technologies, with qPCR generally being less sensitive/more variable.
Question: I am trying to understand how to utilize principal components as covariates for PheWAS studies. Plink can be used to generate eigenvectors, eigenvalues and additionally variant weights. It seems from what I have read that the eigenvectors used as the covariates, but I am unclear then why the eigenvalues are not included in some form. Are the eigenvalues only needed to reconstruct the genetic relationship matrix? And what are the variant weights? Are they the contribution of each variant to each of the principal components?

Answer: Eigenvectors, i.e., principal components, are vectors that are un-correlated and that contain loadings - what they summarise in your data is covariance among your variables. So, when you include eigenvectors as covariates in regression models, you are adjusting your regression coefficients based on these. It has been found that doing this can help to control for population structure in genetic studies. You should decide on the number of eigenvectors (PCs) to include as covariates
via manual inspection of pairwise bi-plots, as per:

- https://www.biostars.org/p/282685/#282691
- https://www.biostars.org/p/271694/

For example, in the following figure, eigenvectors 1, 2, an 3 (PCs 1, 2, and 3) are clearly segregating the different populations; so, we would include these eigenvectors as covariates in modelling if we were studying some disease trait / phenotype that spanned all populations

Kevin
Question: Hi dears,

Can someone tell me a link to plant version of biomart please? I went through this link but it seems inactive http://www.gramene.org/biomart/martview/ea115ee510b63cb6b1680c47ebec3440, in this link http://plants.ensembl.org/biomart/martview/767d21e3edcb8fac82b9d0be5410324b there is not my interest species (citrus). I need plant biomart for id conversion and getting bulk fasta

For my purpose which part I should select from pull-down menu?

http://www.biomart.org/biomart/martview/f6673304b8ab70c4db5ce6b6418f4a96

Answer: Plaza v3.0 is very useful especially in terms of comperative genomics and partly functional annotations.

http://plaza.psb.ugent.be/
Question: I am working on a genome-wide association study. I wish to examine the time to onset of a disease in relation to the alleles present at each SNP (survival analysis). I have done this for a few snps using STATA, but this won't work for a whole genome dataset so I am trying to do it in R.

My files are in PLINK format, and I have included the time to event data in the fam section of the read in plink data using the merge command.

My script is below:

```
#Read in plink files
GWAS<-read.plink("R.bed","R.bim","R.fam")
pheno <- read.delim("C:/Users/USER/Desktop/pheno.txt")
view(pheno)

library("GenABEL", lib.loc="~/R/win-library/3.1")

#find columns
head(GWAS)

#merge plink file with followup and presence of disease
GWAS$fam <-merge(GWAS$fam,pheno, by.x="member", by.y= "Chip")

#this is now included in the plink files
#running in cox model using GenABEL package and mlreg command
coxm <- mlreg(GWAS((fam$followup),(map$allele.1)),"fam$disease")
```

However, it's not working, I'm pretty new at this and I'm not sure how to specify the genotype at each snp for each individual. Has anyone ran a genome-wide cox model in R? If so help would be greatly appreciated!

Answer: I figured it out, thank's everyone for the help. Thought I'd post up the script I used in case anyone else has similar issues.

What my problem was, was that I had the data in the wrong format, so I converted my bed, bim and fam files to ped and map using `--recode` in plink. The files could then be read into R and converted again to the gwaa.data class. Once it was read in the rest was straight forward.

```
#Load package
source("http://bioconductor.org/biocLite.R")
biocLite('GenABEL')

#open package
library('GenABEL')

#set working directory
setwd("C:/where_my_files_are")

#Convert to gwaa.data class
convert.snp.ped(ped='Rel2-0.ped',mapfile='Rel2-0.map', out='rel2-0.out')
survivaldata <- load.gwaa.data(phenofile = "pheno.txt", genofile = "rel2-0.out")

#run quality control
qc1 <- check.marker(survivaldata)
qcGWAS <- survivaldata[qc1$idok, qc1$snpok]

#Run cox proportional hazards model
coxresults <- mlreg(GASurv(survivaldata@phdata$followuptime, survivaldata@phdata$event)~survivaldata@phdata$age, data=qcGWAS)

#View top 20 results
descriptives.scan(results, top=20)
```
Question: Does anyone know if BWA Mem -M flag (mark shorter split hits as secondary) will have a detrimental effect on structural variation calling using split reads with Delly?

I orignially added it to my BWA mapping pipeline for compatibility with Picard/GATK but I believe that issue has been resolved in GATK4. Now wondering if it's even necessary and if it will have an effect on SV calls that use split reads.

Input to this pipeline may include paired, single end, interleaved or split data from illumina, ion torrent, or nanopore.

Answer: For [delly][1], it doesn't make a difference if split reads are marked as secondary or supplementary.

[1]: https://github.com/dellytools/delly
Question: Good afternoon to everybody.
I am working on paired ended whole genome seq data. Specifically I am focusing on structural variant calling and I am at the moment testing different callers.
One of the tools that I am using is Speedseq/Lumpy and, in order to use it, I aligned my fastq files with Speedseq align (which is BWA-based) and generated the three BAM files necessary as input:
- a main BAM file
- a BAM file containing splitted reads
- a BAM file containing discordant reads
All the other callers need a single BAM file as input. So, my strategy was simply to merge into one single BAM file all the three BAM files generated by speedseq. However I found a bit weird the content of the splitted and discordant BAM files, and maybe I am misunderstanding the meaning of splitted and discordant reads.
Here is an example for one read. These are the reads as they come from the fastq files:

R1 fastq file:

`@H0E3UALXX:2:1201:2151855:0 1`

`GTTAGGGTTAGGGTTAGGGTTAGGGTTAGAGATCGGAAGAGCGTCGTGGAGGGAAAGAGTGTTCAGAGCCGTGTAGATCTCGGTGGTCGCCGTATCATTAAAAAAAAAAAAAGGGGGAGTGTGAGGGGGAGGGATCGAGGTAGGGTGCGG`

`AAFFFFJJJJFJJJJJJFJFJJJJJJJFF

Answer: All the alignments are stored in the "main" BAM file, according to the speedseq documentation.

- https://github.com/hall-lab/speedseq#output

The output of SV callers will differ when run on the "main" BAM file and a merged BAM file if the SV caller does not recognize duplicates. By merging the speedseq BAM files, the evidence for any SV is doubled because split and discordant reads are duplicated.
Question: Hi,

I have a question about sharing your code for publication. Do you share all files of your code? I saw some people only share some key functions. Others share the code even with paths to files, but I do not have these files :( . Is there a guide line for sharing code? Thanks!

Answer: There's been a number of articles on the subject, see for example [Ten simple rules for writing and sharing computational analyses in Jupyter Notebooks][1].

[1]: https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1007007
Question: Hi all,

I've been trying to run a sequence alignment using clustal omega from within a python script. However, I keep getting the following error:

FATAL: Sequence file contains no data

However, if I run the raw command on the command line, everything runs fine. Here is the relevant code in my script:

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC
import Bio.SeqIO as SeqIO
from Bio.Align.Applications import ClustalOmegaCommandline

ref_seq = "ERVVIGSKPFNEQYILANMIAILLEENGYKA"
target_seq = "ERVVIGSKPFNEQYILANMINGYKA"

in_file = "unaligned.fasta"
out_file = "aligned.fasta"

# Write my sequences to a fasta file
handle = open(in_file, 'w')
records = (SeqRecord(Seq(seq, IUPAC.protein), id=str(index), name="Test", description="Test") for index,seq in enumerate([ref_seq, target_seq]) )
SeqIO.write(records, handle, "fasta")

clustalomega_cline = ClustalOmegaCommandline(infile=in_file, outfile=out_file, verbose=True, auto=True, force=True)
clustalomega_cline()
#sp.call(str(clustalomega_cline), shell=True) # This will also fail with the same error

**Edit:**

It seems that this is due to either the way biopython writes the FASTA file, or how clustalo processes it.

If I comment out the lines where I write to a file and run the script again, is still fails, even if I run the raw command using subprocess. BUT, if I also make sure that there are no spaces in the FASTA sequence headers, then the commands at the bottom run fine.

I don't know why this isn't a problem when just run from the command line with the headers unmodified. This doesn't make any sense.

Answer: Closing the input file after writing to it, but before calling Clustal Omega, makes this work for me.

Just add

handle.close()

after the call to SeqIO.write(). I think the way it is now, maybe the sequences aren't being written to that file until the script terminates (after the call to Clustal Omega).

Edit: a more elegant version might be:

records = (SeqRecord(Seq(seq, IUPAC.protein), id=str(index), name="Test", description="Test") for index,seq in enumerate([ref_seq, target_seq]) )

with open(in_file, 'w') as handle:
SeqIO.write(records, handle, "fasta")

where the 'with' statement takes care of closing the file handle for you (even if an exception is thrown in the inner block).
Question: Is it possible to load the "broad peak" output files of macs2 ([BED 6+3][1]) into the Integrated Genome Browser (IGB)?

I have tried with no luck, I'm not quite sure whether I need to specify the input format.

[1]: https://genome.ucsc.edu/FAQ/FAQformat.html#format13

Answer: Hi pgarcia,

The .broadPeak format is not a currently supported file format in IGB, however it can easily be converted to the more common .bed format, which is supported.

To convert the file to a bed file, limit the data to the first 6 columns. To do this on the command line use:

$cut -f1-6 FILE.narrowPeak > NEWFILE.bed

Optionally, you can replace the score field (column 5) with any one of the latter fields (columns 7-9 from broad peak format) to make that value appear as the score.

For more information about bed format: https://genome.ucsc.edu/FAQ/FAQformat.html#format1

For more information about broad peak format: http://genome.ucsc.edu/FAQ/FAQformat.html#format12

These file formats (both broadPeak and narrowPeak) have been recently requested by users, and we are preparing to add support for them in the near future.

Nowlan