text:
GCF VS GCA or combination of both for pangenomic studies
I have just started to learn bioinformatics and pangenomics. So if this question seems to you pretty basic then I apologize in advance.
As we know, NCBI for the genome database, there are two kinds of sequences available RefSeq and Genbank. I have read the differences in both RefSeq and Genbank. But I was just curious which one would be preferred for pangenome studies? For example, phylogeny built by RefSeq sequences would be different from GenBank sequences? As I read, RefSeq is already curated, annotated and contaminants are removed. Will it affect the phylogeny that we will build after using any pangenome pipeline?
I am particularly interested in non-synonymous mutations (frameshifts or stop codons inside coding sequences). Do RefSeq curators also remove this kind of unusual errors or mutations? If yes then it would be difficult to study such mutations in RefSeq sequences.
Somebody may say a combination of both Genbank and RefSeq would be good, after removing the duplicate one. Which one should I remove RefSeq or Genbank?.
Regards


text:
Extracting base sequences from ABI/AB1 sanger sequencing chromatogram
I am trying to understand the sanger sequencing ABI/AB1 file format better, and extract base calls from given signal intensities over time.
As I understand, reading in a raw AB1/ABI file into python, I can access the channel corrected values for the different bases in ['DATA9'] to ['DATA12'], and understand which bases belong to a specific channel through ['FWO_1'] i.e either ordered GATC or ATGC etc.
However, the file I am looking at has identified 412 bases as the output from SnapGene, yet has 4950 data records for each channel in the raw ['DATA9'] to ['DATA12']. How do you convert those 4950 data records into the bases correctly? Is there some form of normalisation - e.g every 10 raw records gives data for 1 base position, or peak calling signal detection algorithms. And if so, do you take the average over 10 records, or the highest peak in the 10 records from the channels and the channel with the highest peak/average is the correct base? Do you start at the beginning or do you have some form of cutoff? Does it make sense to convert this figure into a quality score similar to NGS, and do you use the highest peaks in the 10 record interval or the average?
I've tried using scipy's signal.find_peaks() function (python) which can get me various peaks with a minimum height and a minimum distance to other peaks, currently using a minimum peak height of 150 and distance to other peaks of 5, but still getting around 20-30 more peaks than the bases identified.
Could anyone suggest any logic, methods, algorithms or resources that could help with translating the raw signal intensities for the 4 channels into a sequence of accurate base calls?
I know software exists to identify the bases, but that's not what I am looking for. I intend to do this from scratch to understand the workings better so I can write some custom processing algorithms.
I hope this question makes sense.
Thanks

Please consider the answer below because it appears to answer the question and if you are happy to upvote and/or accept it

You can access the peak positions of Sanger Sequencing chromatograms from ["PLOC1"].
The position values in ["PLOC1"] correspond to the index values in ["DATA9"], ["DATA10"], ["DATA11"], ["DATA12"], so you can extract the signal values of each base at the positions.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
Cannot blast database created with blastdb_aliastool: BLAST Database error: No alias or index file found
I want to merge/concatenate BLAST databases as suggested here.
This is my example script that explains what I'm trying to achieve better than words can:
#!/bin/bash

for file in FAM1079.ffn FAM3228.ffn FAM6161.ffn FAM19036.ffn
do
    echo "## Make DB for $file"
 makeblastdb -in $file -dbtype nucl &>/dev/null
    blastn -db $file  -query query_single_nucl.fasta -outfmt 6
    echo
done

echo "## Create testdb.nal..."
blastdb_aliastool -dblist "FAM1079.ffn FAM3228.ffn FAM6161.ffn FAM19036.ffn" -dbtype nucl -title test -out testdb

echo
echo "## Test testdb.nal..."
blastn -db testdb.nal -query query_single_nucl.fasta -outfmt 6

This is the output:
## Make DB for FAM1079.ffn
test_query  FAM1079_00001   100.000 1266    0   0   1   1266    1   1266    0.0 2338

## Make DB for FAM3228.ffn
test_query  FAM3228_000001  100.000 1266    0   0   1   1266    1   1266    0.0 2338

## Make DB for FAM6161.ffn
test_query  FAM6161_000001  100.000 1266    0   0   1   1266    1   1266    0.0 2338

## Make DB for FAM19036.ffn
test_query  FAM19036p_000001    100.000 1266    0   0   1   1266    1   1266    0.0 2338

## Create testdb.nal...
Created nucleotide BLAST (alias) database testdb with 12692 sequences

## Test testdb.nal...
BLAST Database error: No alias or index file found for nucleotide database [testdb.nal] in search path [/path/to/files::]

While the individual databases work, the concatenated databases don't.

The problem was that the database must be specified without suffix, i.e. blastn -db testdb as opposed to blastn -db testdb.nal.
By the way, it produces exactly the same output as
blastn -db "FAM1079.ffn FAM3228.ffn FAM6161.ffn FAM19036.ffn" -query query_single_nucl.fasta -outfmt 6

and the performance is not better.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Problem with SPAdes installation - no bin folder
I have been trying to download SPAdes for an important project I wanted to do, however, I have encountered a small problem while doing so. It appears that there is no 'bin' folder in my SPAdes installation file.
I have used 7-Zip as my main method of unpacking the file (As I am using Windows, I had to download the actual source code rather than a nice executable), and it came in the .tar.gz format, which I believe 7-Zip is rather good with.
A few pictures detailing my problems are shown below.
I hope I can get a little help on this matter. It is slightly time-sensitive, so I would really appreciate a little bit of assistance shortly. Thank you in advance!

As seen above, the first image is the command I used, the usual test command, the second is the file structure I got (Notice the absence of the 'bin' file), and the third is the error I received from running the spades.py --test command.

I see that it supports Linux and MacOS. It appears you are using Windows. Is there a reason you believe that Windows is supported? See https://github.com/ablab/spades/issues/156

No, not at all. I just thought the source code would work on Windows, as there are options for Linux and MacOS downloadables.

Well, sure, 7-Zip will unzip a tar.gz file, but it won't install software for you!

Simply unpacking the source code won't actually install much, you'll need to follow its installation instructions. Note that this might be possible on Windows if you use new bash command line, but in general don't expect any bioinformatics tools to work on windows.


text:
How to convert Bed file to fasta file?
I have Bed file containing start and end of a sequence, and I need to convert it to fasta format, any recommendations?

Please [edit] your question and give us more details. Bed files define regions, do you have a fasta file of the target regions? Is this for a genome? Do you have that genome?

bedtools getfasta can do this. See documentation.

You could start with a folder of samtools faidx-indexed FASTA files for your genome of interest:
$ mkdir /genomes/hg38 && cd /genomes/hg38
$ wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/*.fa.gz
$ for fn in `ls *.fa.gz`; do gunzip $fn; done
$ for fn in `ls *.fa`; do samtools faidx $fn; done

Then use that folder with the following script:
#!/usr/bin/env perl

use strict;
use warnings;
use Getopt::Long;

#
# bed2faidxsta.pl
# --
# Reads BED data from standard input, writes FASTA to standard output.
#
# Dependent on samtools and indexed FASTA data files located in $fastaDir
# variable. Set --fastaDir=dir to set custom directory containing a source
# of per-build, bgzip-compressed FASTA and associated index (fa.gz.fai)
# files, or leave unset to use data in current working directory. Use the
# --fastaIsUncompressed option if the FASTA files are not compressed.
#

# test if samtools is available
`samtools --version` || die "Error: The samtools application is required to run this script. Try 'module add samtools' or install a local copy of samtools.\n";

# default FASTA input is current working directory
my $fastaDir = `pwd`; chomp $fastaDir;
# default is to assume input coordinates use zero-based index scheme
my $oneBased;
# default is to leave IDs alone
my $useIDPrefixAsStrand;
# default is to assume FASTA files are bgzip-compressed
my $fastaIsUncompressed;

GetOptions ('fastaDir=s' => \$fastaDir, 'oneBased' => \$oneBased, 'useIDPrefixAsStrand' => \$useIDPrefixAsStrand, 'fastaIsUncompressed' => \$fastaIsUncompressed);

if (! -d $fastaDir) { die "Error: FASTA directory does not exist\n"; }

while (<STDIN>) {
    chomp;
    my ($chr, $start, $stop, $id, $score, $strand) = split("\t", $_);
if (!defined($chr) || !defined($start) || !defined($stop)) { die "Error: No chromosome name, start or stop position defined\n"; }
    if (!defined($id)) { $id = "."; }
if (!defined($score)) { $score = "."; }
if (!defined($strand)) { $strand = "+"; } else { $strand = substr($strand, 0, 1); }
# adjust coordinates to one-based index, if necessary
my ($queryChr, $queryStart, $queryStop) = ($chr, $start, $stop);
if (!$oneBased) {
        $queryStart++;
}
# adjust strand if required
if ($useIDPrefixAsStrand) {
        $strand = substr($id, 0, 1);
    }
    # lookup
    my $queryFn = "$fastaDir/$chr.fa.gz";
if ($fastaIsUncompressed) {
        $queryFn = "$fastaDir/$chr.fa";
}
my $queryKey = "$queryChr:$queryStart-$queryStop";
my $queryResult = `samtools faidx $queryFn $queryKey`; chomp $queryResult;
# linearize result
my @lines = split("\n", $queryResult);
    my @seqs = @lines[1..(scalar @lines - 1)];
    my $seq = join("", @seqs);
# handle reverse-stranded elements
if ($strand eq "-") {
        $seq = rc_sequence($seq);
    }
    # print to standard output
    my $header = ">".join(":",($chr, $start, $stop, $id, $score, $strand));
print STDOUT $header."\n".$seq."\n";
}

sub rc_sequence {
    my $seq = shift @_;
my $reverse_complement = reverse($seq);
$reverse_complement =~ tr/ACGTacgt/TGCAtgca/;
    return $reverse_complement;
}

Then:
$ perl bed2faidxsta.pl < in.bed > out.fa

This approach can be useful with stranded input, nonstandard genomes, or just to do things with a bit more knowledge of inputs.

Or, since you already have samtools, the poster can use samtools faidx with a file of regions.  It won't be hard to get a bedfile into the right format: http://www.htslib.org/doc/samtools.html

UCSC 'kent' command line tool:
faToTwoBit -bed=yourCoordinates.bed genomeSequence.2bit output.fa

Binary for this command and others available at: UCSC kent command line tools

I think you mean twoBitToFa. Thanks anyway for pointing me in that direction.


text:
How can I identify a recessive and dominant gene?
For example, I have two allelic genes. How can I identify a recessive and dominant gene? Are there any databases with this information? Or the answer to my question is to study the concept of "haplosufficiency"?

I'd probably look up the genes in OMIM, which provides detailed description of the functions and phenotypes associated with genes, including inheritance mode.
For example, this is what you would get for TP53:


text:
What programs account for structural alignment of different parts of distant homologs which have significant structural differences?
If there is a need to perform structural alignment of different parts of distant homologs, which program one should use?
Since distant homologs often have significant structural changes, meaning the assumption of a similar physico-chemical environment for aligned residues would be violated.
Any program/tools or suggested reference will be a great help.
Thanks.

You might want to check out TCoffee's Espresso as it does a structure based alignment (requires structures to be given/exist in PDB). However, in my hands it has always performed less than ideally.

TM-align has some options to do this, see for example https://zhanglab.ccmb.med.umich.edu/TM-align/help.html.
TM-align is similar to TM-score but the alignment is purely structural.

You can use progressiveMauve (http://darlinglab.org/mauve/user-guide/progressivemauve.html), or bubbz (https://www.cell.com/iscience/fulltext/S2589-0042%2820%2930409-0)


text:
Hg38 annotation tracks retrieval
I want to retrieve an annotation object (as GRanges) containing different UCSC annotation tracks such as CpG Islands, TSS, TFBS and promoters for hg38 Human reference genome. But "UCSC Genes" was phased out for hg38. What could be the best resource for this ?
Thanks
RT

The UCSC Data Integrator can combine multiple tracks (CpG Islands, TFBS, cCREs, etc.) into one output annotation file, and it works with any assembly you choose. You can select which output fields you want, then download the tab-delimited file.
From there, you can use R, Python, and/or GNU tools such as awk and sed to coerce the data into a format which is useful for your purposes. You can then load it into R as a data frame and use the columns as arguments in the GRanges constructor to create a GRanges object.


text:
Continuing analysis based on a log-likelihood ratio? (Comparing methylation between two PoolSeq phenotypes)
I'm trying to get some pointers on where to start on the analysis of my dataset. First, some basic info:
Background:
We performed Nanopore sequencing of two poolseq samples called UU and DD referring to damage vs undamaged. We were able to obtain approximately 3x coverage (this is currently more of a proof of concept than a hardcore association job, so the coverage is unfortunate but worth pursuing). Between these samples, we know that one is going to have different methylation patterns in response to a phenotypic trait response, but we don't know what regions along the genome are going to be differently methylated. The goal is to compare what's methylated (as in higher or lower levels of methylation) between the two samples, knowing there ought to be a lot that cancels out and come up with ideally a short list of significantly different regions we can associate to the plastic response.
Output of initial methylation calling data:
Using the raw signal level data against the reference genome and sorted, mapped BAM files the Nanopolish software was used to call methylation. This produces a dataset with the following headers:
1. Chromosome, 2. +/- Strand 3. Start vs. 4. End of the CG dinucleotide (the difference between these is 1 or two bases) 5. The read ID the dinucleotide came from 6. the log likelihood ratio calculated from the model embedded in the methylation calling process using 7. log likelihood that it's methylated and 8. log likelihood that it's unmethylated. 8. number of calling strands (always appears to be 1) 9. number of motifs (don't know about this) and finally 10. the sequence surrounding the dinucleotide, which is always short but variable in length.
I can follow the basic interpretation of the log-likelihood... if it's a positive number, there is statistical support that it's methylated. This would probably be about as much as we need, short of a p-value, if we were looking for a specific region, but we are trying to analyze the differences in de/methylation across the entire methylome. How would you suggest using the statistic values provided in our dataset to compare and visualize the significant different modifications of the dinucleotides between strands/samples?

I assume you don't want to spend 1-2 years developing your own method based on nanopore data, so the path forward is to convert the output into a format used by standard DMR calling programs. There's no standard format for that, annoyingly, but the general things you'll need are:

Coordinates of a C/CpG with coverage
Number of reads supporting methylation at that site.
Number of reads not supporting methylation at that site.

Typically you would do some filtering before/after making such a table so you only have decent quality calls (that would translate into some reasonable log likelihood for your data) and possibly only sites with decent (5x, 10x, 30x, etc.) coverage. You'd then feed that for all of the replicates in all of your groups into programs like bsseq, DSS or Metilene (to not so randomly select 3 examples).
If there aren't DMR finding programs already developed specifically for nanopore data I expect someone is developing it. The technology has a higher error rate, but given that bisulfite conversion as required in Illumina data also leads to a higher error rate that's not necessarily a huge issue if the methylation calling itself is reasonably accurate and the depth is high enough. Note that you probably need 10X coverage as a bare minimum unless your effect size is very very large.


text:
Visual representation of detecting deletion by different read types
I'm having some trouble differentiating between CNV, split reads and read pairs and how each one elucidates the presence of a structural variation, say deletion.
Given a reference and donor sequence, I'd like to see a visual depiction (a link to a reference would be welcome) of a deletion event being detected by each of the above read types. In particular, I'd like to see how each of the read types detects the deletion event itself (i.e structural variant type) and its location (i.e the breakpoint junction).
Bonus points if it also includes some of the other common structural variants like inversion, duplication and insertion.

Welcome to Bioinformatics! This question need more information to be possible to answer. What papers did you read or which tools aren't helpful (and why not?)? Many thanks!


text:
E.coli Sequencing & Analysis
I have been given the task of assembling a 'new' Ecoli genome and analysing the genes present etc.
The Ecoli is a new strain, and has been taken and run on a Nextseq 500 in high-output mode with 150bp paired end reads. The 'raw' files that I have is the forward and reverse reads.
I have initially QC checked the 'raw' files, and subsequently run them through trim galore and checked the QC after that.
For the next step, I now need to assemble my genome. I have been told that SPades will run a 'de novo' assembly for me, and then put that assembly into Prokka for Gene annotation.
Is this the best way to assemble the genome and annotate it? Or should I use another method? I am thinking that I should use a 'mapping' technique to assemble the genome using the Ecoli O157:H7 genome as a reference, but I have no idea how to do this. I would say that I am at an intermediate level with unix, but by no means am I a bioinformatician. Some help and guidance would be greatly appreciated!
Thanks,
Alex.

You can first do assembly of contigs with spades and then use a reference to get longer sequences, for example you can try this: https://github.com/fenderglass/Ragout

A de novo assembly with spades should work as well as a mapping with bowtie2.
While performing a de novo assembly you should also be aware of possible contamination which you may need to remove. By using a mapping this may not be necessary.
When you use a mapping, you must be aware that you may miss relevant and diverged regions because of the mapping to a reference strain. So a mappping would be preferrable if your target genome is closely related to your reference genome.
I'd like to give you the advice to use some kind of distribution system like Anaconda which helps you to keep your system clean and make it easier to reproduce your steps.
Here are some possible steps for both ways. I am using an e.coli k12 sample, sra number SRR6940089 :
De novo assembly
1. Quality Check ( you already did this )
2. Run Spades
spades.py -o k12 -m 8 -t 4 --pe1-1 SRR6940089_1.fastq --pe1-2 SRR6940089_2.fastq &> k12.log

This command will start an assembly using 4 threads and 8 GB RAM. You can scale it to whatever your system is capable of or what you'd like.
You now have a folder called k12 which holds the assembly information.
To get some information about the contigs you can use the emboss toolkit.
infoseq k12/contigs.fasta

fasta::contigs.fasta:NODE_859_length_78_cov_1.000000 -              NODE_859_length_78_cov_1.000000 -              N    78     43.59
fasta::contigs.fasta:NODE_860_length_78_cov_1.000000 -              NODE_860_length_78_cov_1.000000 -              N    78     44.87
fasta::contigs.fasta:NODE_861_length_78_cov_1.000000 -              NODE_861_length_78_cov_1.000000 -              N    78     43.59

This is the end of the contigs.fasta. So you see contigs.fasta contains 861 contigs. And you also see thate ther length and coverage is quiet low. So you need to remove this low coverage regions.
3. Remove low quality contigs
There are multiple ways how you could achieve this. I'd like to filter sequence smaller than 500 bp. the tricky part is to identify the coverage threshold. Contigs with a low coverage could be some artifacts or contamination.
4. Check for possible contamination
So to check for possible contamination you can use kraken2. But it needs some disk space. I am using here the minikraken database
kraken2 --db minikraken --threads 4 --output ecoli.kraken --report ecoli_kraken.report k12/contigs.fasta &> kraken.log

The ecoli_kraken.report shows an overview of how the contigs could be classified to availible species from the database.
Based on this results you can also remove some contigs.
5. Run Prokka
Now after you removed unwanted contigs, you can run prokka to get the annotation.
prokka --cpus 4 --outdir prokka k12/contigs.fasta &> prokka.log

And you will find the information in the prokka folder
Mapping
The mapping process is a little bit easier. You can use bowtie2 and samtools
bowtie2 -x ecoli -1 SRR6940089_1.fastq -2 SRR6940089_2.fastq -p 4 | samtools view -@ 4 -O BAM -o ecoli.bam

158216 reads; of these:
  158216 (100.00%) were paired; of these:
    35696 (22.56%) aligned concordantly 0 times
    120047 (75.88%) aligned concordantly exactly 1 time
    2473 (1.56%) aligned concordantly >1 times
    ----
    35696 pairs aligned concordantly 0 times; of these:
      4458 (12.49%) aligned discordantly 1 time
    ----
    31238 pairs aligned 0 times concordantly or discordantly; of these:
      62476 mates make up the pairs; of these:
        56512 (90.45%) aligned 0 times
        5365 (8.59%) aligned exactly 1 time
        599 (0.96%) aligned >1 times
82.14% overall alignment rate

You now also get some information about the similarity to your reference sequence.
Now you need to work a little bit more on the BAM file. You need to sort and index this.
samtools sort ecoli.bam -o ecoli_sort.bam
samtools index ecoli_sort.bam

Now you can create a consensus sequence using bcftools and seqtk
bcftools mpileup -f sequence.fasta ecoli_sorted.bam | bcftools call -c | vcfutils.pl vcf2fq | seqtk seq -aQ64 -q 20 -n N > ecoli_consensus.fasta

This gives you now a sequence where unknown bases are marked as N.
As I wrote before, using anaconda would make the installation a lot easier. So you can write the following lines into the file bio.yaml
name: bio
channels:
  - bioconda
  - defaults
dependencies:
  - aragorn=1.2.38=h470a237_2
  - barrnap=0.9=2
  - bcftools=1.9=h4da6232_0
  - bedtools=2.27.1=he941832_2
  - blast=2.6.0=boost1.60_0
  - bowtie2=2.3.4.3=py36h2d50403_0
  - emboss=6.6.0=h6debe1e_0
  - hmmer=3.2.1=hfc679d8_0
  - infernal=1.1.2=h470a237_1
  - kraken2=2.0.7_beta=pl526h2d50403_0
  - libdeflate=1.0=h470a237_0
  - libidn=7.45.0=2
  - minced=0.3.2=0
  - parallel=20160622=1
  - perl-app-cpanminus=1.7044=pl526_1
  - perl-bioperl=1.6.924=4
  - perl-carp=1.38=pl526_1
  - perl-constant=1.33=pl526_1
  - perl-data-dumper=2.161=pl526_2
  - perl-encode=2.88=pl526_1
  - perl-exporter=5.72=pl526_1
  - perl-extutils-makemaker=7.34=pl526_2
  - perl-file-path=2.15=pl526_0
  - perl-file-temp=0.2304=pl526_2
  - perl-parent=0.236=pl526_1
  - perl-threaded=5.22.0=13
  - perl-xml-namespacesupport=1.12=pl526_0
  - perl-xml-parser=2.44=pl526h3a4f0e9_6
  - perl-xml-sax=1.00=pl526_0
  - perl-xml-sax-base=1.09=pl526_0
  - perl-xml-sax-expat=0.51=pl526_2
  - perl-xml-simple=2.25=pl526_0
  - perl-yaml=1.27=pl526_0
  - prodigal=2.6.3=1
  - prokka=1.12=3
  - samtools=1.9=h8ee4bcc_1
  - spades=3.12.0=1
  - tbl2asn=25.3=0
  - boost=1.60.0=py36_0
  - bzip2=1.0.6=h14c3975_5
  - ca-certificates=2018.03.07=0
  - certifi=2018.10.15=py36_0
  - curl=7.62.0=hbc83047_0
  - expat=2.2.6=he6710b0_0
  - fontconfig=2.13.0=h9420a91_0
  - freetype=2.9.1=h8a8886c_1
  - giflib=5.1.4=h14c3975_1
  - icu=58.2=h9c2bf20_1
  - jpeg=9b=h024ee3a_2
  - libcurl=7.62.0=h20c2e04_0
  - libedit=3.1.20170329=h6b74fdf_2
  - libffi=3.2.1=hd88cf55_4
  - libgcc=7.2.0=h69d50b8_2
  - libgcc-ng=8.2.0=hdf63c60_1
  - libgd=2.2.5=hceca4fd_3
  - libpng=1.6.35=hbc83047_0
  - libssh2=1.8.0=h1ba5d50_4
  - libstdcxx-ng=8.2.0=hdf63c60_1
  - libtiff=4.0.9=he85c1e1_2
  - libuuid=1.0.3=h1bed415_2
  - libwebp=1.0.0=h222930b_1
  - libxml2=2.9.8=h26e45fe_1
  - ncurses=6.1=hf484d3e_0
  - openjdk=8.0.152=h46b5887_1
  - openssl=1.1.1=h7b6447c_0
  - perl=5.26.2=h14c3975_0
  - pip=18.1=py36_0
  - python=3.6.7=h0371630_0
  - readline=7.0=h7b6447c_5
  - setuptools=40.5.0=py36_0
  - sqlite=3.25.2=h7b6447c_0
  - tk=8.6.8=hbc83047_0
  - wheel=0.32.2=py36_0
  - xz=5.2.4=h14c3975_4
  - zlib=1.2.11=ha838bed_2

Now you can run following command and most of the installation part is completed. Just the database for kraken is missing.
conda env create -f bio.yaml

I hope this gives you some help how to process your data. If you need some information just let me know ;)


text:
What could cause differing counts of R1 and R2 in Paired End Sequencing (RNASEQ)
I recently finished mapping an RNAseq run using STAR2.3.0 and noticed the read1 and read2 counts are different, according to samtools flagstat.
The map% is ~80% but the R1 R2 counts are:
R1=139283692
R2=137472495
This usually does not happen for PE sequencing, which makes me a bit nervous to proceed with analysis. What are some potential reasons something like this could occur? I do not have experience with this happening and can't seem to find reasons by searching the web.

Presumably, the actual read counts in your two fastq source files are the same, right? I have no experience with RNAseq but the first thing that comes to mind is a problem in the sequencing step where too many reads were of low quality and were discarded by the aligner. Does that make sense?

The actual read counts of the two fastq source files are the same. It makes sense and was a thought of mine. I was wondering if there is a concrete way to check this out?

Use FastQC to compare the overall quality of both fastq files.

@d_kennetz Does the actual read counts matches R1? Have you checked which of the reads are missing in the bam file?


text:
Extracting all reads from bam file which match read IDs in another file
I have a long list of read IDs of interest to me in a file called read_names.txt. it is simply in the format:
m54197_180831_211346/4981510/ccs
m54197_180831_211346/6226723/ccs
m54197_180831_211346/6619607/ccs
...

etc where these are the actual read names from a fastq file. I am then trying to find those in a bam file for which the fastq has been mapped, and I can accomplish this like:
for ID in `cat read_names.txt`
do
samtools view inbam.bam | grep $ID >> read_locs.sam
done

However, this method is obscenely slow because it is rerunning samtools view for every ID iteration (several hours now for 600 read IDs), and I was hoping to do this for several read_names.txt files.
I tried sort of flipping the script a bit and running samtools view first but it only returned the first read ID present in the file and stopped:
samtools view inbam.bam | for ID in `cat read_names.txt`; do grep $ID >> read_locs.txt; done

Am I missing something or is this the only way to accomplish this task?

Just as a general rule never do for var in $(cat file). Also known as bash pitfall number 1.

thanks terdon, I will keep this in mind and find the better ways to handle this in the future. I use this method quite often so that is great for me to know!

Yeah, most of us use it very often until it bites us! :)

It is still slow but grep has a -f option to take in a file
samtools view inbam.bam | grep -f read_names.txt > read_locs.txt

This took 30 seconds as opposed to several hours, so I'd say the overhead was definitely in samtools view. I feel silly because of how simple this was. Thank you!

just as a warning this can take much longer if the read name file is big.

I appreciate the heads up! Comparatively, I do not think it will ever take longer than the method I was using so it is still a great answer. (This is serving as a caveat to your warning).

samtools can do this natively too using -N, --qname-file CLI option:
samtools view -N read_names.txt in.bam > read_locs.sam

❯ samtools --version
samtools 1.16.1
Using htslib 1.16


text:
Albacore basecalling running but outputs 0 reads
I am trying to basecall data produced by the MinION using the SQK-LSK109 kit and FLO-MIN106 flowcell via the command-line. My version of albacore is the latest (v2.3.4). I tried running using the following command-line (the name of the sequencing run and fast5s are different and are default from machine output):
python3.6 /home/dkennetz/.local/bin/read_fast5_basecaller.py -i ./seqrun/fast5/ -t 1 -s ./good_reads/ -f FLO-MIN106 -k SQK-LSK109 -r -o fastq

-t 1 is just for this test. It is running and finishing almost instantly with 0 errors but outputting no reads. I have previously run this same command-line using the SQK-LSK108 kit and had no problems.
I then went on to validate the fast5s using the h5_validate tool:
python3.6 /home/dkennetz/.local/bin/h5_validate /home/dkennetz/.local/lib/python3.6/site-packages/h5_validator/schemas/single_read_fast5.yaml seqrun/fast5/0/seqrun.fast5 -v

The output from this is:
    Matcher <'UniqueGlobalKey', type=group match_type=exact count=
{'minimum_count': 1, 'maximum_count': 1}> was not satisfied after matching <HDF5 file "seqrun.fast5" (mode r)>

Error at /:

Matcher <'PreviousReadInfo', type=group match_type=exact count=
{'minimum_count': 1, 'maximum_count': 1}> was not satisfied after matching <HDF5 file "seqrun.fast5" (mode r)>

I read on the community page that this may be because the header is missing some relevant meta-data information? Any help on how to fix this would be greatly appreciated!
This is cross-posted on the community page for ONT here:
https://community.nanoporetech.com/posts/albacore-basecalling-runni

I believe you have multiread fast5 files, as generated by recent versions of MinKNOW. A conversion script is available to convert these to the older, one-read one-file format. AFAIK Albacore is being deprecated and no version will be made supporting the multiread fast5 files, as Albacore is being replaced by Guppy.
The command line to convert from multi read fast5s to single read fast5s is as follows:
$multi_to_single_fast5 -i /path/to/fast5/ -s /path/to/save/output/

Hi Wouter, thanks for the answer! This worked as expected. To clear any ambiguity, I edited to add the command-line for future viewers. You can choose to accept or not, I just thought it added clarity! I appreciate the link.


text:
Is there a read mapper that allows input in SAM format instead of FASTA/FASTQ?
When we get raw sequenced reads in FASTQ format, we usually need to do some pre-processing (QC, demultiplexing, etc.) prior to sequence alignment. The only way to register the results of this preprocessing is to add metadata to each read. This is best done by using the SAM format which has a well-defined system of metadata tags.
However, it seems that not a single read mapper (except maybe for bwa aln with option -b) accepts SAM (or one of the binary equivalents BAM/CRAM) as input to make use of the metadata that has been added during preprocessing. This implies that one has to convert the unaligned reads back to FASTQ format.
In any case, whether one uses FASTQ throughout or converts back to FASTQ, a common annoyance with FASTQ is that metadata is registered as part of the query name - this is not standardized or even well-defined. After alignment, we again have to do some text processing to get the metadata out of the query name and into proper metadata tags defined by SAM. This back and forth hacking is highly unfortunate and error prone. I have seen this, for example, in some single-cell sequencing projects where we have highly multiplexed sequencing runs with reads that need to be QC'd and classified. The new SAM standard from May 2018 does indeed define proper tags for handling cell-multiplexed runs. SAM is the de-facto standard to store both unaligned and aligned reads for years to come so it should be standard for read mappers to accept it as input.
I would be grateful for some hints and if we could collect here some read mappers that do accept SAM/BAM/CRAM as input format.

Why not store the read identifiers and metadata in a separate table file?

Of course, this is possible. But apart from not being an answer to my question, this is not a good solution because it is not standardized and invites inconsistencies. You have to write special code to read the table. The SAM format (which is also a table) is maybe not the ideal but it is the best we have. If we could all stick to the SAM format instead of introducing custom formats the interoperability of bioinformatics tools would be much better.

I’m a bit confused because the preprocessing you mention is pretty routine, and is generally done without the involvement of BAM files.

Please remember to confirm an answer once you've received one.

With bwa-mem or minimap2, the recommended way is
samtools fastq -T BC,RX name-grouped.sam | bwa mem -C -p -

This passes the BC and RX tags through bwa-mem and copies them to the SAM output. If you are running paired-end alignment and your BAM is coordinate sorted, you have to run samtools collate to group reads by name:
samtools collate -Ol0 in.bam tmp | samtools fastq -T BC,RX - | bwa mem -pC -

Here bwa mem -C or minimap2 -y copies FASTQ comments to SAM output. This strategy is much easier to deal with than SAM/BAM. A caveat is that it doesn't work well for read groups because each read group requires a separate header line. You have to manually modify the header.
I am not sure which short-read RNA-seq mappers have implemented a similar feature. You might consider to send a feature request.

The common solution for scRNA-seq is to put cell barcodes and such in read headers and then post-process things with UMItools.
But regarding your actual question, STAR can accept SAM/BAM as input with the --readFilesType SAM PE option (it's SAM for both that format and BAM). You can swap SE for PE if your data is single-end (or is effectively that way, as is the case in many scRNA-seq experiments). Note that in the case of PE experiments that mates must be next to each other. I really don't know what this does with auxiliary tags though, so you would need to see if they're lost before you start processing a lot of data).

At the Broad Institute, we do something similar to what is described by user172818. We use uBAM (unmapped BAM) to store unmapped reads, because as you say the SAM format has a more appropriate way to encode metadata. However, we use a Picard tool called SamToFastq to do the initial conversion then pipe the output of that to BWA mem. Then we use an additional step to reattach metadata, read group info and so on, using a Picard tool called MergeBamAlignments.
You can see the full implementation here (starts as line 71):
https://github.com/gatk-workflows/gatk4-genome-processing-pipeline/blob/master/tasks/Alignment.wdl

SNAP is a hash-index based genomic aligner that can take SAM and BAM as input. It has been used as an aligner for RNA-seq (see SNAPR) and BS-seq (see BS-Seeker3). SNAPR seems to accept SAM/BAM but it is not clear whether BS-Seeker3 wraps this functionality as well.
I tried SNAP and it seems to do the job well, i.e. it only adds tags during alignment and does not remove metadata that has been deposited prior to alignment, neither from the SAM header nor from individual reads.
NOTE I post this answer to my own question after I became aware of SNAP.


text:
Can someone recommend software for reference-based genotype-estimation?
Problem: I have low-coverage sequencing data from humans in vcf format. For each sample, I have PL scores / preliminary genotype calls from the software ATLAS (it's a bit like GATK for those who don't know). Because the samples are low coverage (mean 1x, with lots of missing data), I am looking to enhance the genotype calls by using a reference panel and genetic map, before go on to phase the samples. Briefly, the reference panel and genetic map can be used as a kind of prior to inform one about genotypes when there is little data present.
What I've done so far: So far I have been using Beagle4 with a command:
beagle4 gl=lowCoverageSamples.vcf ref=highCoverageReference.vcf map=geneticMap out=lowCoverageSamples.newCalls.vcf

However, Beagle4 is extremely slow, very memory inefficient and quite old. There have been many advances in phasing methods, which are now much more accurate, faster and memory efficient than the beagle4 model. However, I haven't come across any similar improved methods for the genotype re-estimation step.
Does anyone knew of any newer pieces of software which do such a thing, which may more faster and more accurate than Beagle4. It needs to have 2 things in particular:
1) Accept data input in PL score format, as to allow the appropriate representation of uncertainty in low coverage samples.
2) Produce a genotype call / genotype probability for each SNP in each individual as output.

In case anyone was interested, a new piece of software for rapidly imputing and phasing low coverage sequence data from genotype likelihoods has just been released. It is called GLIMSPE and is exactly what I was after.
It looks to be much faster (1200x speedup) compared to Beagle 4.0, and is more accurate at imputing low frequency variants. It is also more amenable to use with larger reference panels (such as the Haplotype Reference Consortium panel).


text:
What are fundamentals to work on bioinformatic data analysis ?
Well, I'm a MSc student in Information Technology Engineering, trying to become more familiar in the Complex networks specially Bioinformatic field to get more knowledge of cancer detection by Biochemical data analysis.
beside python / R language what else should I work on to prepare myself for an academic seminar?

You need to be much more specific in order to get helpful answers. What analysis are you planning to do and which tools are you currently using? If you want a broader overview of the field, there are many courses, books, and blogs on this topic.

Which academic seminar (About which topic? what are the ponents?) ? What do you know about biology and health? Are you familiar with experiments performed on the laboratory? BTW this question is quite broad, we can answer about a concrete answerable questions. You could try asking at biostars.org but please add more information to the question in any case


text:
How to create a graph of multi column gene dataset?
well, I'm new to Bio. what I'm going to do is create a network/graph of a dataset related to the genes I think. here's the sample of the nodes dataset with CSV format named lung-nodes :

and this is a sample of edges I think with trrust name :

these data are stored in two separate file
what I must do is create a graph(network) of these nodes and edges by Python and a visualization library like networkX, but I do not how.
I must make a relation between these to the dataset by creating a graph of them
can you help me out?

You must work in Python? Is this a homework assignment?

@benn, yes it is. I am not familiar so much with gene and protein properties. I already create a graph from social network dataset but they weren't like this. they had just two column in a single file. but right now i got two separate CSV file which contains bio info which i must create a graph of them. i think R studio could make it easier

As far as I understand, it seems that the first table you provided contains biological/experimental information related to your genes, i.e. the nodes of the network you will build. The second table is the one you need to build the network: each line is an edge. You can use it as input and build the network using the IGraph library which works for both Python (as you requested) or R (as you suggested as an alternative).
Most probably, but I can't say for sure, the first table was just a dataset from which you (or someone for you) recovered a list of nodes that where then submitted to a database, i.e. String, commonly used to recover information related to the interactions between your nodes, i.e. genes. Once you have this info (in the second table) you can build the network!
edit:
if you're not familiar with IGraph, you can have a look at this question describing how to obtain an IGraph graph from a csv file or this other one which uses a different csv file format to build a graph


text:
Finding the conserved region of protein in a given set of Kingdom
We have to find the conserved region of proteins in a given set of kingdom and we have to give more weightage to distantly related organism. What approach I should take to solve this problem?

Welcome to the site. In order to be able to answer the question better, could you explain what data/information do you have? Did you search for some terms related to your question? Why what you found didn't help? What else have you tried?


text:
What are 2D reads in the Oxford MinIon?
Reading the Mash: fast genome and metagenome distance estimation paper, I see the reference to "2d reads" notion under the minION context (a sequencing technique)
e.g
    In both cases Mash was able to
correctly differentiate these closely related species
(ANI ≈ 95 %) using 43,806 and 91,379 sequences collected
from single MinION R7.3 runs of B. anthracis
Ames and B. cereus ATCC 10987, respectively (combined
1D and 2D reads) <----------  HERE.

But couldn't find any reference in Wikipedia or online papers to what it isCan someone shed some light?

see here: https://www.biostars.org/p/176781/

Early MinION sequencing runs had forward and reverse DNA templates joined together by a hairpin adapter, so that the sequencer would read both strands from the same template. The consensus sequence that was generated by combining the base calls of opposite strands is referred to as a 2D read.
Due to a dispute with Pacific Biosciences, and for other reasons related to chemistry and basecalling, the equivalent mode in current sequencing runs is called 1D². There is no physical ligation of forward and reverse templates, but the kit and flow cell chemistry are modified to encourage the reverse strand to follow the forward strand very soon afterwards through the pore, with the joining and consensus generation carried out in software after the sequencing is done. In the output fast5 files, these 1D² reads can also be referred to as 2D.


text:
Purpose of correlating docking score with catalytic efficiency(Kcat/km)

What is the purpose of correlating docking score with catalytic efficiency (kcat/km)?
I saw some papers are doing this, is it to validate the docking protocol?
Lets say the docking protocol is correct, the docking score will positively correlate with kcat/km?

Usually the Km I've seen related to K_cat is the michaelis-mendel constant, not a docking score (but maybe this is different terminology). Could you explain some context of this question? A text of the paper or a reference where do you find this. Also if you could show what you have searched/read it would allow us to skip these resources when trying to find the answer.

https://www.nature.com/articles/ncomms12965

refer to the graph in page 3

What is the purpose of correlating docking score with catalytic efficiency (kcat/km)?

This is specific to the enzymology of paper you are citing. If the solubility of your substrate is lower than your Michaelis constant (K_M) you cannot fit to a Michaelis-Menten curve and can only get the slope of the initial straight part of the curve, which is the catalytic efficiency (k_cat/K_M). Ideally K_M would be better. But in generally catalytic efficiency ought to correlate with well calculated ddGs of ligand binding.

I saw some papers are doing this, is it to validate the docking protocol?

Totally. Docking finds were ligands bind. Although the energetics don't match too well the experimental kinetics as of 2019.
The catalytic efficiency (k_cat/K_M) is the ratio of the turnover (k_cat) over Michaelis constant (K_M). The Michealis constant is in turn the ratio of on-rate over off-rate plus turnover. The ΔΔGibbs of the ligand ought to correlate with the ratio of the on-rate/off-rate, while the turnover is unrelated.
The turnover should correlate with the difference in energy between the substrate bound form and the transition state. This is the realm of theozyme design and requires correctly protonated ligands and protein, but when it is you get a nice transition state hill (e.g. personal example).
There are a few "ought to" in the above. Docking and force-field energy scores are not precise. The paper you refer to is doing very advance docking (akin to dynamic undocking) and, not accidentally, the ligands are very flexible and rather soluble.

Lets say the docking protocol is correct, the docking score will positively correlate with kcat/km?

A positive correlation coefficient just means that generally both the catalytic efficiency and docking scores increases or decrease together. If a metric "worked backwards", such as ΔΔGibbs (more negative is better), you would want a negative correlation. What you want is a strong correlation, with the sign a result of the nature of your metrics.


text:
Post docking evaluation
What to do after protein-ligand docking? Is there any necessary or enhancement step recommended after docking to validate the confirmation of highest score is the best? For your information, I am using CCDC GOLD for docking and GOLD score as my scoring function.

Why do you need to to a docking? Do you have a biological question? For example, I want to know if this drug docks to this receptor. To validate the docking you usually test it in real life...

If you really have to, you can try different docking algorithms and scoring FF. However, you might get subtle differences due to how your ligand is parameterised.

just a follow up to the answer, what if we don't want to proceed with the assays and just want to stick with the bioinformatics. What else can I do to make my filtered drug molecule a stronger candidate?

You can run a Molecular Dynamics simulation to see if the pose is stable and whether the protein changes/adapts it’s conformation.
You can calculate the ligand efficiency, the docking score divided by the number if heavy atoms, and rank the compounds accordingly.
You can also compare the scores of known ligands to the experimental values and see if the docking the scores are predicting known ligands or correlate with the experimental ones.
At some point you have to make a selection and purchase or synthesize the compounds. How many compounds can you test experimentally? How many can you purchase? What is the purpose of these ligands going to be? Drugs? Did you screen the correct set of ligands e.g. a library of lead compounds that are ready to purchase for cheap, that would be a good start.

To further this answer, a form of steered MD that I have found promising is (dynamic undocking)[https://www.nature.com/articles/nchem.2660]: your ligand gets yanked out. The output is not super clear though and requires some manual checking.


text:
How to find Wild Type (WT) protein structure on Protein Data Bank (PDB)?
Currently, I am working on Galactose oxidase enzyme for protein docking and I am totally new to it. I have to study how mutation on specific residue will make mutated galactose oxidase better than the wild type. My questions:

There are so many structure available on PDB, how to know which structure is the so called wild type structure? Is the structure with PDB Code: 1GOG belongs to wild type?
I was requested to extract the ligand(substrate) out of the structure and docked it back. But seem all the structures available on PDB do not complexed with ligand. For what I know, the substrate for galactose oxidase could be D-galactose. Could I use the galactose oxidase structure from pubmed and dock it to the enzyme? Btw, what is purpose of extract the ligand out and dock it back?

Welcome to the site. 1) Did you check out the description of the structures? Of which organism are you interested in the wild type of  the protein? 2) Yes you could. Why do you want to extract the ligand from a structure and dock it back? (This seems homework and could be only to you to practise how to use the software)

Yes, I did read the literature and it mentioned that 1GOF is belong to wild type but interestingly, I found another 2 structures (1GOG, 1GOH) that have the identical sequence as wild type (1GOF). That's why I confused which structure to use, btw the differences between these 3 structures are only their ligands. My question is since their sequences are identical, does it means that I could just use either one of them to perform docking?

All three structures are from the same specie but other species could have a different wild type. There is no rule why a single molecule can't be linked to several others but the structure might be different (or change when docked with other ligands). Yes you could use any structure to perform docking but the suggested structures could be different according to the initial state of the structure


text:
How to interpret the GROMACS MD Simulation RMSD Graph?
Hello brilliant researchers!
After virtual screening, I performed MD Simulation via GROMACS to check the stability of the protein-top hit complex by following this tutorial http://www.mdtutorials.com/gmx/complex/09_analysis.html
After that, I facing some difficulties in interpreting the resulted RMSD graph as it is quite different from most of the published RMSD graph.
Based on the graph attached, how to interpret it? Is AChE-dexketoprofen complex (Red) more stable (or conserved) than BuChE-dexketoprofen complex (Black)?

The opposite-ish.
Your graph is showing only the ligand, but RMSD is generally used for structures as a whole —so I'll mention about whole protein for future readers.
RMSD is a measure of how much does a structure deviate in position from a reference. Your ligand in the protein in red "wiggles" more than the one in black.
Structure RSMD
If this were RMSD of the protein (and not the index group for the ligand), it would not really mean that the black protein is more stable. Stability is to do with force-field energetics not position. Some protein have two state positions, say MBP is a flexing dumbbell. This would be seen as a jump in the RSMD graph, similar to that seen for the black trace.
However, the command you used was based on gmx rms -s em.tpr -f md_0_10_center.xtc -n index.ndx -tu ns -o rmsd_jz4.xvg from the page you linked. Key here is the index -n index.ndx, which limits it to the ligand. Before this happened gmx trjconv aligned the protein in the trajectory, so that common mistake was not done.
Ligand RSMD
In terms of ligand RMSD, lower deviation/less wobble correlates with more stable binding, but energetics (gmx energy) will tell you how strongly does it bind.
The black trace has a jump, so your ligand in the black trace (which wobbles less) has an alternative binding mode. This is why you are looking at RSMD.
To explore this further you can:

look at the trajectory
in PyMOL using the command split_states twice, once for the first and then the last state, followed by align of the two
If it were the whole protein, looking at the per bond matrix -bm option in Gromacs rms and spotting the change.

3Å??
I strongly suggest checking the trajectory because both have a deviation of 3Å, which for a ligand is not good.

Did the energy minimisation converge properly?
If it is because the proteins themselves wobble substantially, revisit gmx trjconv with a focus on the active site minus the ligand.
Is the ligand binding but leaving one substituent group to wobble?

Energetics
As with docking, the "score" is free energy difference (in kcal/mol). This is not the same as k_off and k_on and their ratio K_d (Dissociation constant). For those you need to do steered MD or have a way longer simulation (not recommended). The reason being is that the ligand has to make its way into the active site. However, you have a secondary binding mode in black, so have two free energy scores, this means that there is a stronger k_on.


text:
Error in formation of correlation matrix-(corrupt matrix -- dims not not match length)
I am working on a script to create a correlation matrix and mutual_rank matrix from RNAseq gene counts for 57992 genes in 7027 samples. I have already tested this script for a smaller file (for 10,000 genes) and it is working, however, when I am running this script on my final input file, it is giving the following error:
 corrupt matrix -- dims not not match length (this error comes after gene to
 gene correlation step)

Script which i am running is:
E <- read.table("mdat_human_test.txt", header = TRUE, row.names=1)

# gene-to-gene correlation
C <- cor (t(E))
write.table(C, "human_correlation_genes_mdat.txt", quote=FALSE, sep="\t",
row.names=TRUE)

# correlation rank matrix
R <- t(apply(1-C, 1, rank))
write.table(R, "human_rank_genes_mdat.txt", quote=FALSE, sep="\t",
row.names=TRUE)

# mutual rank matrix
M <- sqrt(R * t(R))
write.table(M, "human_mutual_rank_genes_mdat.txt", quote=FALSE, sep="\t",
row.names=TRUE)

q()

What is your biological question of interest ? Why are you doing this correlations? In addition to Devon's answer, I would like to point out that if you are using all the human annotated genes you might need to filter and normalize them. You'll presumably will have 25% of "genes" which are poorly described and/or with few counts are you really interested in keeping those ?

You can't compute a correlation matrix when you have more than ~46000 rows, since a standard R matrix can have a maximum of 2^31-1 values. Have a look at packages like bigcor. Alternatively, consider if having 3+ billion correlation values is really useful.


text:
How to fit all genes (labels) in chromosome ideogram plot made by RCircos package?
I am using Rcircos to make a chromosome ideogram plot for my gene list (n = 45). However, I am getting this error:
 Not all labels will be plotted.
 Type RCircos.Get.Gene.Name.Plot.Parameters()
 to see the number of labels for each chromosome.

I am using the below-given script and what to know is there any way to fit all my genes in the chromosome ideogram plot
 library(RCircos)
 data(UCSC.HG19.Human.CytoBandIdeogram)
 RCircos.Set.Core.Components(
 cyto.info=UCSC.HG19.Human.CytoBandIdeogram,
 chr.exclude = NULL,
 tracks.inside = 10, tracks.outside = 0)

 #Making a plot with RCircos
 out.file <- "RCircos_IF.pdf";
 pdf(file=out.file, height=8, width=8, compress=TRUE);
 RCircos.Set.Plot.Area();
 par(mai=c(0.25, 0.25, 0.25, 0.25));
 plot.new();
 plot.window(c(-2.5,2.5), c(-2.5, 2.5));

 # plot chromosome ideogram
 RCircos.Set.Plot.Area();
 RCircos.Chromosome.Ideogram.Plot()
 gene.lable.data = read.table ("IF.txt", sep="\t", header=T, row.names = 1) # IF.txt is my input file
 RCircos.Gene.Connector.Plot(gene.lable.data, track.num = 1, side = "in");
 track.num <- 2;
 RCircos.Gene.Name.Plot(gene.lable.data,name.col = 4, track.num)

 dev.off()

My input file looks like
    Chromosome  chromStart  chromEnd    Gene.name
    1   chr1    150363091   150476566   RPRD2
    2   chr1    150549369   150560937   ADAMTSL4
    3   chr1    91949371    92014426    BRDT
    4   chr1    31365625    31376850    FABP3
    5   chr1    150960583   150975004   CERS2

To change the maximum number of genes in RCircos, we could modified the char.width according to this link.
The actual value in a RCircos session could be modified with get and reset methods for plot parameters.
params <- RCircos.Get.Plot.Parameters() #$char.width
params$char.width <- 100 #default 500
RCircos.Reset.Plot.Parameters(params)

#the maxLabels are updated accordingly
RCircos.Get.Gene.Name.Plot.Parameters()

Just want to add that apart from params$char.width, we can also use
params$text.size
params$text.size <- 0.2 #default 0.4
params$base.per.unit < 50000 #imp #default value = 30000


text:
How to fix the Error in Kruskal Wallis rank sum test for comparing the guides in python? Also to check whether using the right test for my data?
df
gene_name   guide_1 guide_1_new Correlation
MMP-1   A   A   1
MMP-1   A   B   0.426115
MMP-1   A   C   0.522499
MMP-1   A   D   0.431587
MMP-1   B   A   0.426115
MMP-1   B   B   1
MMP-1   B   C   0.60113
MMP-1   B   D   0.534858
MMP-1   C   A   0.522499
MMP-1   C   B   0.60113
MMP-1   C   C   1
MMP-1   C   D   0.622206
MMP-1   D   A   0.431587
MMP-1   D   B   0.534858
MMP-1   D   C   0.622206
MMP-1   D   D   1

I used four guides in an experiment (data = df). I calculated correlation between them based on the their log fold changes. I want to check which guide is outperforming or underperforming in comparison to other three guides in my experiment. I tried Kruskal and Wallis test to do this analysis with the script given below. However, I got an error ValueError: Samples must be one-dimensional. Also is this the right test to do this analysis?
My script
from scipy import stats

x = df[['guide_1']]
y = df[['guide_1_new']]
z = df[['Correlation']]
stats.kruskal(x, y, z)

I am a little confused, you are passing the columns of guide_1 and guide_1_new as samples to scipy.stats.kruskal? Those are non-numeric data. And the correlations, while numeric, suggest to me that you are doing something a little different from what Kruskal-Wallis is implementing. It looks like you might be trying to pass in the categories as vectors of strings, but I don't think it works that way (?). And I don't think applying it to correlation coefficients will be very interpretable. Not an expert on scipy.stats but this doesn't look right. I'd suggest going back to raw data.

Answered below, if you are stuck I can talk you through the calculation if you perform a df.to_dict() and post the output.

Reformatting the data for the KW and doing the test is easy.
Using "guide_1" alone and pivot() to get A = [1, 0.426115 ...
df2 = df.pivot(columns='guide_1', values='Correlation')

This will now present A, B, C, D etc ... into columns...
x = df2[['A']].dropna()
y = df2[['B']].dropna()
z = df2[['C']].dropna()
a = df2[['D']].dropna()

stats.kruskal(x, y, z, a)

That will work ... whether that is a good test is a complicated question because this is log transformed data. Transformation is used in parametric tests. Removing the log transformation is strongly recommended because KW works with raw data.
You can switch the calculation to guide_1_new instead. It will give a different answer ... but I don't know the experiment. The original guide looks cool.
Looking at the data - the test should work because I believe KW is 'unpaired' (you need to check) and looking at the value the distribution looks pretty similar between A, B, C and D. However, just because the KW says 'yes' doesn't mean thats a robust answer (log data). If the untransformed data says 'yes' ... thats a robust answer for KW. Just check the 'unpaired' stuff - thats important for the result.

Personal note, I originally used a third df but was not needed,
df3 = pd.concat([df2[col].sort_values().reset_index(drop=True) for col in df2], axis=1, ignore_index=True)

I'm a little confused here, it looks to me in the OP like the numeric input is correlation coefficients from comparing A to A, A to B, etc. So for A to A, you get an unsurprising correlation of 1, which is now included in the A group, which doesn't seem like a very informative input. Does OP really care about whether guide_1 group A is more correlated to other groups than those groups are among themselves (which is what this is testing for, effectively)? Doesn't make sense to me, shouldn't OP go back to raw (pre-correlation) data? Probably missing something here.

@MaximilianPress KW will give a result representing the data. Agreed non-parametrics prefer raw data. The presence of A->A , B->B will of course bias the test, but the gist is A ~ B ~ C ~ D - and will not be 'thrown out of court'. Capturing log-change or non-log change isn't easy. A 2-way ANOVA type (4 columns vs 2 rows) would do this - absolutely, but its complicated and non-parametrics are not good at this calculation (Friedman test, statisticians avoid it). Using 1-way ANOVA on this specific data would be a mistake. Getting it 'spot on' is an investment. Right now 'it's ok'.


text:
Transcription factor targets for Human and Mouse genes
I want to get transcription factor targets for my list of genes (For Human and Mouse)
I have found ENCODE Transcription Factor Targets database (The link is given below). But I am looking for some other databases also if available. I will really appreciate if anyone can suggest any other database.
http://amp.pharm.mssm.edu/Harmonizome/dataset/ENCODE+Transcription+Factor+Targets
Thank you

What is the problem with ENCODE? I am aware of the existence of TRANSFAC, but I think is a non-free resource

@Plat My only problem with ENCODE is that I am not sure when was it last updated (I think in 2012)

JASPAR is another useful transcription factor profiling database.

The JASPAR CORE contains a curated, non-redundant set of profiles, derived from published and experimentally defined transcription factor binding sites for eukaryotes. It should be used, when seeking models for specific factors or structural classes, or if experimental evidence is paramount.

http://jaspar.genereg.net/

JASPAR is only based on TF motifs, which are very short sequences and not enough to precisely tell you the true TF-motif relationships, especially when multiple TFs have complex logic gates regulating overlapped targets.
An alternative way is to find some transcriptomic data upon perturbation of each TF. The drawbacks are that the redundancy of TFs will mask a lot of targets, and that the up/down regulated genes may not be direct targets but targets of targets (of targets of targets...).
Another way is to find ChIP-seq data for TFs you care about, and hopefully you can find the right biological context. But still ChIP signal isn't exactly binding pattern, which in turn isn't exactly indicating target relationship.
Overall, it's a challenging task to know what are the true targets of a TF.


text:
Map gene IDs to protein structure if exist
I want to map some gene ids to their proteins structures if they exist.
The genes are rows of this CSV. There are other information that may help also such as Product description but that may lead to anbiguities. I am looking for a programmatic way, biopython or pdb API?

From EMBL SIFTS you can download a table that goes from Uniprot ID to PDB and chain and numbers that can give you any offset. Do note that offset correction is a problem for about 10% —I use python's pymol2 module to fix these (yup, silly but best option). Also if the organism is common you might wanna check out SwissModel's premade threaded models of close homologues(20% id min)

Not sure if you want the structure exactly or a handle of some sort, but the best way i know to get ALL your 5000 identifiers to map to uniprot/pdb is uniprot's mapping service.
You can either copy your ids straight from excel or parse with python ex.:
import pandas as pd
import os

file_name = 'aap7847_Zhang_SM_Table-S5.xlsx'
sheetname = 'Table S5'
df = pd.read_excel(file_name, sheet_name=sheetname, skiprows=1)
ids = list(df['Gene_ID'])

name_of_file_to_save_to = "geneIds.txt"
with open(name_of_file_to_save_to, "w") as newf:
    for id in ids:
        print("{}\n".format(id), file=newf)
    newf.close()

You can then save the Uniprot's response and either add it as a column to your file or map to other parts of Elixir ecosystem with SIFTS, for example, as @Matteo Ferla says.

In the case that you maybe need a single id's information, you could hit their API too, although this wouldn't work for bulk or at least didn't work for me. You can map the uniprot accession you get from this to PDB or whatever else you need!
import requests
def get_uniprot_accession(geneid, format='list'):
    resp = requests.get(
        "https://www.uniprot.org/uniprot/?query={}&sort=score&format={}".format(geneid, format)).text
    if format == 'list':
        resp = resp.strip()
    return resp


text:
failed to install hdf5r
I failed to install hdf5r using this command provided by https://satijalab.org/seurat/mca_loom.html: devtools::install_github(repo = 'hhoeflin/hdf5r'),below is the code,can you help me ?


text:
Separate boxplots for multiple violin plot
I am using the following function from seurat package to generate multiple violon plots and I am interested in adding box plots to them but it doesn't work when I have plotted different data at once.
Is there a way to solve it ?
For example, this works:
library(Seurat)
VlnPlot(object = pbmc_small, features.plot = 'PC1') +
  geom_boxplot()

But this will simply lead into an empty box on top of my plots:
VlnPlot(object = pbmc_small, features.plot = c('PC1', 'PC2')) +
  geom_boxplot()

This happens because the violin plots are combined using cowplot::plot_grid before being returned by VlnPlot. You can prevent the plots from being combined by setting combine=FALSE, then modify each one by adding a boxplot, then combine the modified plots using Seurat::CombinePlots.
However, the combine argument is currently broken in VlnPlot. I have pushed a fix but it's not on the public branch yet. VlnPlot is just a wrapper around ExIPlot (expression by identity plot) in Seurat v3 so right now this will work:
library(Seurat)
library(ggplot2)

plots <- ExIPlot(object = pbmc_small,
                 features = c('PC1', 'PC2'),
                 pt.size = 1,
                 combine = FALSE)

for(i in 1:length(plots)) {
  plots[[i]] <- plots[[i]] + geom_boxplot() + theme(legend.position = 'none')
}
CombinePlots(plots)

Note this answer is for Seurat v3

The vioplot package comes built in with boxplots. You can download it from CRAN or there are more features (including formula input and separate colours) in the development version on GitHub: https://github.com/TomKellyGenetics/vioplot
devtools::install_github(“TomKellyGenetics/vioplot”)
library(“vioplot”)
vioplot(pbmc_small@dr@pca[,1]~pbmc_small@meta.data$ident)

To plot several PCs:
library(“RColorBrewer”)
colours <- brewer.pal(length(unique(pbmc_small@meta.data$ident)), “Set3”)
par(mfrow = c(1, 2))
vioplot(pbmc_small@dr@pca[,1]~pbmc_small@meta.data$ident, main = “PC1”, xlab = “Identity”, col = colours)
vioplot(pbmc_small@dr@pca[,2]~pbmc_small@meta.data$ident , main = “PC2”, xlab = “Identity”, col = colours)
par(mfrow = c(1, 1))

You can add points to these with jitter: plot(pbmc_small@dr@pca[,1]~jitter(pbmc_small@meta.data$ident, 0.5), add = TRUE)

This version is now support on CRAN


text:
amplicon sequencing set up
I have .fastq files from a targeted amplicon sequencing by Illumina MiSeq run. I want to get the target sequences to perform Blast search on them. What is the pipeline to use for processing amplicon sequencing data? I also have 20% PhiX.

Usually you have the targeted regions based on the library/ panel design prior to sequencing.  Do you not have that information?  It is standard Illumina processing to remove the PhiX reads before creating the FASTQs

This is crossposted on biostars: https://www.biostars.org/p/357824/

Do you want the reads that contain the amplicon? If you know the primer sequences you could do this pretty easily


text:
VEP output SIFT_score unclear
We have been experimenting with VEP (Variant Effect Predictor). One of the meta data attributes that we are interested in is the SIFT score, indeed when we apply the dbNSFP plug we get a column containing the scores (named SIFT_score). However, I don't understand why there are sometimes dots or multiple values in the fields. For example, the gene ENSG00000196924 below has 5 transcripts:

The SIFT_score column contains several values, not 1 per transcript/rs-number...
Here is another example that confuses me (I added the SIFT_pred column this time):

There are two mutations, the lower one can be expressed in 4 transcripts and I thus understand that there can be 4 SIFT scores, but why are all for in given in every row? Is the first one the SIFT_score for the first transcript?
One last example, again 4 transcripts, but now 2 of the scores are dots, what does that mean?

I have been looking for quite some time now how to interpret this data, any help is appreciated.

The first gene you mention, ENSG00000196924, actually has 6 transcripts (link to the VarSome.com page of variant rs371839875), not 5. It's just that one of them is non-coding:

So the Sift scores you see are indeed one per transcript, it's just that there are 6 because dbNSFP also includes a score for the non-coding transcript of the gene.
The dots are just there as placeholders, they mean there was no value associated with that transcript. Many tools will show some sort of symbol instead of an empty field both for clarity and for practical technical reasons.
Visiting the variant's page on VarSome gives you a clearer picture since we collapse the identical scores and also include the converted rankscore provided by dbNSFP so you can have a single number for your variant:

Disclaimer: I work for the company behind VarSome, but it's a free tool. You need to pay to annotate VCF files (unlike the 100% free VEP), but it's free to use as a lookup tool for single variants.

The dbNSFP plugin from VEP accesses tables of data for each variant from dbNSFP and pulls out the values. dbNSFP provide their SIFT scores in that format: a score for every transcript affected by the variant, all on one line. The lookup is just for the variant, not for the variant/transcript combo, so they provide scores for all variant/transcript combos. You can also get a column that gives you a list of the transcripts or proteins (Ensembl_transcriptid or Ensembl_proteinid) in order so you know which score goes with which transcript.
A better way to get SIFT scores with VEP is to get them directly from VEP, rather than using dbNSFP. This will get the SIFT score that goes with the transcript on the line with the relevant transcript.

I am guessing the dots are there for cases where dbNSFP doesn't have a value for the relevant transcript, right?

Yes, that's it. Could be that the variant isn't missense in that transcript.


text:
Snakemake: Migrating from deprecated cluster.json to new profiles.yaml
I am an avid user of Snakemake. Recently we have been refreshing our pipelines and I saw that a cluster.json file is no longer the recommended way to store the cluster configuration.
I used to start my pipelines like this:
snakemake --cluster-config cluster.json --cluster "qsub -W umask=002 -d cluster_outputs  -l nodes={cluster.nodes}:ppn={cluster.ppn} -N {rule}" --jobs 245 --use-conda --rerun-incomplete --keep-going --latency-wait 60

With the cluster.json looking like this:
{
    "__default__" :
    {
        "cluster"   : "qsub",
        "jobs"      : 70 ,
        "nodes"     : 1,
        "ppn"       : 1
    },
    "rsem_quantify" :
    {
        "ppn"       : 8
    },
    "star_map_human" :
    {
        "ppn"       : 4
    },
}

Now, I understand that some of this goes in to the new profiles.yaml like this:
cluster: qsub
jobs: 70

But how do I set the ppn values per rule? And how do I set the other configuration options into profile.yaml (like: snakemake "qsub -W umask=002 -d cluster_outputs  -l nodes={cluster.nodes}:ppn={cluster.ppn} -N {rule}", etc.)
Any examples are appreciated, https://github.com/snakemake-profiles is not really helping

Just wanted to note that I asked the same question in stackoverflow. No  good answer yet though.

It looks like you can just bundle your cluster config in its own YAML file, and reference that from within a profile YAML.  Here's a clumsy little example using a fake job submission script: https://github.com/ressy/example-snakemake-profiles

The current documentation makes it sound like cluster configuration support is going away and being replaced by profiles:

While still being possible, cluster configuration has been deprecated by the introduction of Profiles.

...but from those example profiles, that might be misleading; here's an example template for a profile YAML that references a separate cluster-config YAML:
https://github.com/Snakemake-Profiles/generic/blob/master/%7B%7Bcookiecutter.profile_name%7D%7D/config.yaml
...
cluster-config: "$((INSTALDIR))/overwrite_cluster_param.yaml" #abs path
...

That cluster-config YAML has the same default and per-rule entries just your JSON version does:
https://github.com/Snakemake-Profiles/generic/blob/master/%7B%7Bcookiecutter.profile_name%7D%7D/overwrite_cluster_param.yaml
...
# specific parameters for certain rules, which need more time/memory

#run_assembler:
#  queue: bigmem
#  time: 100
#   threads and memory definen in config file

So I think you can just set up a pair of YAML files accordingly including options for your rsem_quantify and star_map_human rules.  From my own test it looks like it handles the cluster argument templating just as it did before when using a profile.  (If this is correct we should nag them to make the new documentation more explicit though.)


text:
Get SRR numbers from a GSE number
I'm trying to programmatically download files from SRA (with prefetch in sra-tools) starting from a GSE number, using Python.
There is an explanation for using entrez-direct on the CLI, using esearch and the efetch. But I want to use Python. This is what I have now:
from Bio import Entrez
Entrez.email = 'my_name@email.me'
handle = Entrez.esearch(db='gds', query='GSE134195', term="biopython")
record = Entrez.read(handle)

which gives:
record = Entrez.read(handle)
record

{'Count': '17', 'RetMax': '17', 'RetStart': '0', 'IdList': ['200066260', '200045602', '200037064', '200030751', '100013939', '300762998', '300762997', '300762996', '300762995', '300762994', '300762993', '300762992', '300762991', '300762990', '300762989', '300762988', '300762987'], 'TranslationSet': [], 'TranslationStack': [{'Term': 'biopython[All Fields]', 'Field': 'All Fields', 'Count': '17', 'Explode': 'N'}, 'GROUP'], 'QueryTranslation': 'biopython[All Fields]'}

But what are these numbers? IdList looks useful, but there are 17 numbers where I expect 14. How do I get to SRR (SRA) files? The files I want can be found manually here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE134195, click on SRA run selector (at the bottom), this bring you here: https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA554267, there I can download the Accession list that I want. But I'd like to automate that process.
Thank you in advance.

I think thats easy, if I get chance I'll put together a bit of code.

That would be very nice of you, thank you.


text:
Python string editing in Snakemake
I'm running a custom perl script using Snakemake. I use this rule:
rule complexity_20mer_counter:
    input:
       os.path.join(fastq_trimmed_dir, '{sample}_{read}' + fastq_extension.split('.')[0] + '_val_' + '{read}'.strip('R') + '.fq.gz')
    output:
        os.path.join(custom_qc_dir, '{sample}_{read}_20mer_counter.txt')
    shell:
        '''
        #!/bin/bash
        perl scripts/20mer_counter.pl {input} > {output}
        '''

It fails (I'll get to why later) with this error:
Building DAG of jobs...
MissingInputException in line 1 of /rst1/2017-0205_illuminaseq/scratch/swo-406/2019-01-24/rules/20mer_counter.smk:
Missing input files for rule complexity_20mer_counter:
fastq_trimmed/xs0067_P2018SEQE27R01_S1_R1_001_val_R1.fq.gz

... this is because I use the output of trim_galore, this rule looks like this:
rule trim_galore:
    input:
       [os.path.join(fastq_raw_dir, '{sample}_'+ read + fastq_extension)
       for read in reads]
    output:
        timmed_fastq_files = [os.path.join(fastq_trimmed_dir, '{sample}_' + read + fastq_extension.split('.')[0] + '_val_' + read.strip('R') + '.fq.gz')
        for read in reads],
        trimming_reports = [os.path.join(trim_galore_qc_dir, '{sample}_'+ read + fastq_extension + '_trimming_report.txt')
        for read in reads]
    conda:
        "../envs/fastqc.yaml"
    shell:
        '''
        #!/bin/bash
        trim_galore --fastqc --gzip -o {fastq_trimmed_dir} --paired --retain_unpaired {input}
        # Move all qc reports to the fastq_trimmed_dir
        mv {fastq_trimmed_dir}/*.zip \
           {fastq_trimmed_dir}/*.html \
           {fastq_trimmed_dir}/*.txt \
           {trim_galore_qc_dir}
        '''

As you can see, trim_galore not only uses my {read} parameter in the output filename but again adds it at the end, only now without the "R".
I tried solving it in my first rule by applying a python .strip('R') on my {read} wildcard, but this does not seem to work. How can I correctly get '1' and '2' from 'R1' and 'R2' (which are my {read} wildcards)?
(Extra bonus points for a nice output.trim_galore.something as input for my complexity_20mer_counter rule, which does not work because the output is a list (it needs 2 fastq files at the same time), but my complexity_20mer_counter rule can be run in parallel again.)

For me it's quite hard to give you a final solution, without having a look on the whole Snakefile or at least a runnable minimal example. I would suggest that the values for read are already without the R. Than you can easily decide when a R is needed in the filename and when not.

This Snakefile is quite difficult to follow, and its complexity means that subtle mistakes might easily creep in. For example, in your trim_galore rule you have timmed_fastq_files instead of trimmed_fastq_files as the label for your output files. I don't know if this is a problem, but little things like that can definitely mess up the whole pipeline. Also, the shebang at the beginning of each "shell" directive is unnecessary. :-)

Yes, I already caught the timmed error, it was not it. I did try putting reads to ['1', '2'] and making _R part of the string but may people use my pipelines and I want to make the default what every is used most. Doing what @finswimmer suggests will lead to a lot of code to make the "R" optional somehow. Although one could probably figure something out. I would just like to do string operations in the wildcard, like {read}.strip('R'). I'll try to look into making something minimal working but that is a lot of effort in this large pipeline.

It seems that it is not possible to apply strip() directly on the wildcard placeholders. Luckily it is possible to executed functions in the input directive.
So try changing your input of the complexity_20mer_counter rule to this:
lambda wc: os.path.join(fastq_trimmed_dir, wc.sample + '_' + wc.read + fastq_extension.split('.')[0] + '_val_' + wc.read.strip('R') + '.fq.gz')

The wildcard object is passed to our anonymous function. So have access to the wildcards value via wc.sample and wc.read, which are now handled as normal strings and we can apply strip on them.

Yes, this is what I was looking for, it works, thank you! (I did call "wildcards", not "wc", not sure if that is supported, I call "wildcards" everywhere.


text:
Running Snakemake in one single conda env
I have been experimenting a lot lately with Snakemake, I love it. Recently I also switched to using conda (--use-conda) in the way that is advertised. However, I have some issues with it, mostly related to the way we work. I work in a lab that processes data from our sequencer but also external datasets, moreover, multiple people use my pipeline.
So, with respect to conda this means: currently I have 5 envs, averaging 1.5 GB per env, this takes about 20 min to set up and it's quite some data to pull down. Especially since it is pulled down every time someone analyses data (sometime multiple times a day). This feels like quite a waste of resources. Also, I fix some package versions in the yaml files, but certainly not all (this may lead to issues when repos get reshuffled I read), so I also can't guarantee the exact environment long term (say over years, dependencies may change/update).
So I'm considering switching to a central conda env which I version and of which I export lists of all packages installed for reproducibility reasons. This will then contain all tools.
My question is, is there a nice way to do this? I had issues with putting a "source activate" line in every shell command (and how to do this for Python code under "run" anyway) before, see this: https://stackoverflow.com/questions/54024862/snakemake-conda-env-parameter-is-not-taken-from-config-yaml-file.
So, is there a way, also when using a qsub cluster, to always use a specific conda env, for everything you do?

Using only the --use-conda directive, snakemake copies environments somewhere in the local .snakemake directory. This is repeated for each directory you run it in and replicates all environments effectively for each dataset.
You can however tell snakemake to use a shared environment location using --conda-prefix as described here:

Specify a directory in which the ‘conda’ and ‘conda-archive’
  directories are created. These are used to store conda environments
  and their archives, respectively. If not supplied, the value is set to
  the ‘.snakemake’ directory relative to the invocation directory. If
  supplied, the –use-conda flag must also be set. The value may be given
  as a relative path, which will be extrapolated to the invocation
  directory, or as an absolute path.

This way you could only download/copy the environment once and share it, e.g. this way:
snakemake --use-conda --conda-prefix /path/to/shared/environment

Yes, this is what I was looking for, great, thanx. In this way I can still create multiple environment but they are created once and available to everybody.

--conda should be replaced by --use-conda. Former is not supported by snakemake as far as I could tell.

Your are absolutely right, this has probably changed in new versions. I've replaced it and also the link was updated. Thanks for the pointer

A good practice is to provide yml file along with your Snakefile. This file should contain all the enviroment definitions that are needed to run your pipeline. The minimum content of this file would be to install snakemake itself. Any other program used in your pipeline should be included as well. For reproducibility version numbers must be used.
So, if someone would like to use your pipeline, (s)he first creates a new enviroment based on your definition:
$ conda env create -f my_pipeline.yml

Afterwards activate this environment and run the pipeline:
$ conda activate my_pipeline
$ snakemake

The conda directive in your Snakefile is only neccessary if you need to use program which dependencies cannot be resolved within the main environment. If you start snakemake with the --use-conda parameter, snakemake will create a new environment based on the provided definition file and activate it for the rule.

Especially since it is pulled down every time someone analyses data (sometime multiple times a day).

The environment creation is only done once. If you run the pipeline a second time snakemake can activate the environment diretcly without downloading and installing something again.
The tutorial I wrote some time ago on biostars about snakemake and conda might be useful for you.

Ah, sounds good, I see that "run:" parts also don't support a conda env, so you really need the main env to have everything as well. I'm testing this now, thanx. Will mark as accepted answer when everything works.

Inside your Snakefile you can provide the same conda: directive to each rule. If you install your pipeline in a central location and have everyone use the same conda when running it (hint: use a wrapper script that sets the path to conda) then everyone will be using the same environment. This has the added benefit that you can make this readable to everyone, but disallow them to make any modifications (so users can't mess things up).
The tack we've taken is to have a basic conda env with snakemake and little else and then separate envs for various pipelines (and versions thereof if software versions in your YAML files change) or parts thereof, in cases where you have rules using tools with conflicting dependencies. This ends up being fairly convenient for more complicated pipelines, though if you can get away with just a single environment then finswimmer's answer is simpler.
Note that all of this should work with any scheduler, since snakemake will insert the environment activation stuff in the script it submits to your scheduler.
As an aside, this is exactly how snakePipes works.


text:
Snakemake: Re-read the samples.tsv file midway through the workflow
All, I have a workflow that in some rules requires the sample's read-lengths. For this I always ran a script that would parse all fastq files and would add the read lengths to the sample table. Rules could then take the read lengths for each sample from there during execution.
However, I recently made a workflow/pipeline that downloads files from SRA, so I need to compute the read lengths at a later stage during workflow execution (after prefetch/fasterq-dump).
I added a step to update the sample table, however, I think the sample table is only read at the very beginning of the workflow.
What would be the nicest way, the Snakemake preferred way, to re-read the samples.tsv file before a rule is executed?
In general I like storing meta data in the sample table so, I'll probably use this mechanism more often in the future, if I find a nice methods to re-parse the samples table.
Right now I do this in rules:
rule deletion_profile:
    input:
        rseqc_input_bam,
        rules.add_max_read_length.output.done_file
    output:
        os.path.join(rseqc_dir, '{sample}.deletion_profile.txt')
    conda:
        "../envs/rseqc.yaml"
    params:
        read_length = lambda wildcards: samples.loc[wildcards.sample, 'read_length']
    shell:
        '''
        deletion_profile.py -i {input} -l {params.read_length} -o {rseqc_dir}/{wildcards.sample}

Thanx.

You're right, it doesn't re-run things again for each job, so your original definition of samples is still the one it's working with.
It seems like it goes against the grain of Snakemake to overwrite a file in-place that acts as input to a rule, so I think the most Snakemake-y way you could do it would be to treat your samples.tsv as part of the input/output graph like everything else.  You could have a rule with samples.tsv as input and samples.somethingelse.tsv as output where that latter one has read lengths and whatever else gets automatically added.  (Alternatively, you could just have your updated samples.tsv read within the rule itself rather than relying on a shared copy loaded earlier, which is kind of halfway to the above idea.)
Either way, maybe something like:
rule deletion_profile:
    input:
        bam=rseqc_input_bam,
        samples="samples.extended.tsv"
    output:
        os.path.join(rseqc_dir, '{sample}.deletion_profile.txt')
    conda:
        "../envs/rseqc.yaml"
    run:
        samples = whatever_loads_samples(input.samples)
        read_length = samples.loc[wildcards.sample, 'read_length']
        shell("deletion_profile.py -i {input.bam} -l {read_length} -o {output}")

(On a side note, you can use {output} rather than repeating yourself for that part.)


text:
Trying to use GEO provided RNA-seq counts but provided Human gene annotation table is incomplete
I recently learned that GEO has been processing RNA-seq raw read into counts (for Differential Expression analyses) and TPM values using a HISAT2 based pipeline. This is attractive because one can save a lot of time and compute using these values. The only issue is that I sometimes need to compare the values to local data (processed with RSEM and more recently with Salmon) or to data that has not been processed by GEO (yet?).
So I started comparing GEO's data with my own. To give some background: I have always been using Ensembl genomes and annotations but with switching to Salmon I started using Gencode as Salmon devs recommend it.
To compare I processed the set GSE223024 myself and I compared the results with the GEO provided values (see here).
When investigating the TPM values I was not able to directly compare the values because the GEO data uses NCBI Gene IDs (ints), whereas I always used Ensembl IDs (ENSG... numbers). There is an Human gene annotation table (Human.GRCh38.p13.annot.tsv.gz), it contains columns with (NCBI) GeneIDs, Hugo Symbols and EnsemblGeneIDs, however only about 3/4th of the genes (27077 out of 39376) have an Ensembl ID in this file. For example NCBI ID 653635 (WASH7P) does not, but looking online it does have an Ensembl ID (ENSG00000227232).
This leaves me with a problem, I need to translate all NCBI GeneIDs to Ensembl IDs but the provided file is incomplete.
I have some questions:

How to best translate ALL NCBI geneIDs into Ensembl IDs?
Are these IDs always 1 to 1 translatable into Ensembl IDs?
Why is this annotation file incomplete?
Is the GEO pipeline consistent? I mean is it always the same pipeline with the same references?

How do you all use the GEO values in your pipelines?

Please don't put multiple questions into the same post. If you have different questions, submit them as separate posts. In particular, your questions about the GEO pipeline seem separate/distinct from your questions about gene database conversion.

If you want to compare with a particular gene database, your best approach will likely be to map directly to that database. Different gene databases will always have different ideas about what is considered to be a gene (otherwise they wouldn't be different databases).
As mentioned a few times elsewhere, translating between different databases is never going to be a 1:1 mapping, so should be avoided if remapping is possible:
Converting Gene Symbols to Entrez ID in ambiguous cases

There will be some gene IDs that have no associated Ensembl IDs
There will be some Ensembl IDs that have no associated gene IDs
There will be some gene IDs that have multiple associated Ensembl IDs
There will be some Ensembl IDs that have multiple associated gene IDs

If you can't do that with your own data, or don't want to do that, then you should re-map the GEO data from the original raw FASTQ files.