text:
Efficiently iterate over bed and subset pyranges objects
I am using pyranges and I have two datasets that I want to compare in several genome-wide intervals based on a bed file.
For every interval in the bed file, I want to get all overlapping positions from the two datasets and sum two columns. That gives me 4 numbers for a fisher exact test, but that's not the most important part.
My current code looks like below, and just iterates over the itertuples() of the bed file as pandas dataframe, and uses that to slice the pyranges object:

import pandas as pd
import pyranges as pr
from pyranges import PyRanges
from scipy.stats import fisher_exact
import numpy as np

def main():
    colnames = ["Chromosome", "Start", "End", "num_motifs", "calls", "methylated", "freq", "seq"]
    h1 = PyRanges(pd.read_csv("dataset1.tsv.gz", sep="\t", names=colnames, header=0))
    h2 = PyRanges(pd.read_csv("dataset2.tsv.gz", sep="\t", names=colnames, header=0))
    for chr, begin, end in pr.read_bed("chr21.bed.gz").merge().as_df().itertuples(index=False, name=None):
        odr, pv = fisher_exact([meth_counts(h1[chr, begin:end]),
                                meth_counts(h2[chr, begin:end])])

def meth_counts(interval):
    if not interval.empty:
        all, methylated = interval.as_df()[['calls', "methylated"]].sum()
        return [all - methylated, methylated]
    else:
        return [0, 0]

This code works well, but I'm not sure if I use the "full power" of what pyranges could offer me for iterating and subsetting. I uploaded some example data on google drive (250GB).

Note that gr[chr, begin:end] will create a new NCLS to do lookups for every single iteration. This is sloooow. It is better to use functions like overlap, intersect, join etc.

Setup:
import pandas as pd
import pyranges as pr
from pyranges import PyRanges
from scipy.stats import fisher_exact
import numpy as np

# ! zcat dataset1.tsv.gz | head -2
# chromosome    start   end num_motifs_in_group called_sites    called_sites_methylated methylated_frequency    group_sequence
# chr21 5010053 5010053 1   3   0   0.000   CACCACGTCCA

# avoid using more columns than we need, especially seq which contains object
# is heavy to lug around
colnames = ["Chromosome", "Start", "End", "calls", "methylated"]
usecols = [0, 1, 2, 4, 5]

h1 = PyRanges(pd.read_csv("dataset1.tsv.gz", sep="\t", names=colnames, header=0, usecols=usecols))
h2 = PyRanges(pd.read_csv("dataset2.tsv.gz", sep="\t", names=colnames, header=0, usecols=usecols))

regions = pr.read_bed("chr21.bed.gz")
regions = regions[~regions.Chromosome.str.contains("_")].merge()

regions.ID = list(range(len(regions)))

regions
# +--------------+-----------+-----------+-----------+
# | Chromosome   | Start     | End       | ID        |
# | (category)   | (int32)   | (int32)   | (int64)   |
# |--------------+-----------+-----------+-----------|
# | chr21        | 5065660   | 5065810   | 0         |
# | chr21        | 5068140   | 5068290   | 1         |
# | chr21        | 5128400   | 5128550   | 2         |
# | chr21        | 5154780   | 5154930   | 3         |
# | ...          | ...       | ...       | ...       |
# | chr21        | 46637040  | 46637190  | 2377      |
# | chr21        | 46638480  | 46638630  | 2378      |
# | chr21        | 46639820  | 46639970  | 2379      |
# | chr21        | 46692680  | 46692830  | 2380      |
# +--------------+-----------+-----------+-----------+
# Unstranded PyRanges object has 2,381 rows and 4 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

h1
# +--------------+-----------+-----------+-----------+--------------+
# | Chromosome   | Start     | End       | calls     | methylated   |
# | (category)   | (int32)   | (int32)   | (int64)   | (int64)      |
# |--------------+-----------+-----------+-----------+--------------|
# | chr21        | 5010053   | 5010053   | 3         | 0            |
# | chr21        | 5010215   | 5010215   | 5         | 1            |
# | chr21        | 5010331   | 5010335   | 6         | 0            |
# | chr21        | 5010525   | 5010525   | 2         | 2            |
# | ...          | ...       | ...       | ...       | ...          |
# | chr21        | 46699713  | 46699718  | 274       | 254          |
# | chr21        | 46699739  | 46699739  | 133       | 113          |
# | chr21        | 46699756  | 46699756  | 37        | 32           |
# | chr21        | 46699770  | 46699773  | 156       | 130          |
# +--------------+-----------+-----------+-----------+--------------+
# Unstranded PyRanges object has 331,373 rows and 5 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

Get total number of methylated per region:
def count_total_and_methylated(df):
    """For each region in bed, get sum of calls and methylated"""

    grpby = df.groupby("ID")
    total = grpby.calls.sum()
    methylated = grpby.methylated.sum()

    chromosome = df.Chromosome.iloc[0]
    start = grpby.Start.first()
    end = grpby.End.first()
    region_id = grpby.ID.first()

    return pd.DataFrame({"Chromosome": chromosome, "Start": start, "End": end, "ID": region_id,
                         "calls": total, "methylated": methylated})

rh1 = regions.join(h1).drop(like="_b").apply(count_total_and_methylated)
rh2 = regions.join(h2).drop(like="_b").apply(count_total_and_methylated)

rh1.File = 1
rh2.File = 2

def add_missing(df, df2):

    """Since join is inner rh1 and rh2 might not have the exact same intervals. Fix this."""

    missing_in_first = np.setdiff1d(df2.ID, df.ID)
    first_missing = df2[df2.ID.isin(missing_in_first)].copy()

    f = df.File.head().iloc[0]

    first_missing.loc[:, "File"] = f
    first_missing.loc[:, "calls"] = 0
    first_missing.loc[:, "methylated"] = 0

    return pd.concat([df, first_missing])

rh1 = rh1.apply_pair(rh2, add_missing)
rh2 = rh2.apply_pair(rh1, add_missing)

j = rh1.join(rh2, suffix=2).drop(like="File*|ID*|(Start|End)2")
j
# +--------------+-----------+-----------+-----------+--------------+-----------+---------------+
# | Chromosome   | Start     | End       | calls     | methylated   | calls2    | methylated2   |
# | (object)     | (int32)   | (int32)   | (int64)   | (int64)      | (int64)   | (int64)       |
# |--------------+-----------+-----------+-----------+--------------+-----------+---------------|
# | chr21        | 5065660   | 5065810   | 188       | 8            | 70        | 0             |
# | chr21        | 5068140   | 5068290   | 168       | 166          | 74        | 72            |
# | chr21        | 5128400   | 5128550   | 901       | 14           | 143       | 1             |
# | chr21        | 5154780   | 5154930   | 653       | 0            | 593       | 0             |
# | ...          | ...       | ...       | ...       | ...          | ...       | ...           |
# | chr21        | 46637040  | 46637190  | 66        | 63           | 69        | 62            |
# | chr21        | 46638480  | 46638630  | 28        | 0            | 30        | 3             |
# | chr21        | 46639820  | 46639970  | 22        | 0            | 18        | 14            |
# | chr21        | 46692680  | 46692830  | 71        | 25           | 57        | 23            |
# +--------------+-----------+-----------+-----------+--------------+-----------+---------------+
# Unstranded PyRanges object has 1,922 rows and 7 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

Compute fisher exact using fisher instead of scipy:
def fe(df):
    from fisher import pvalue_npy

    cols = "calls methylated calls2 methylated2".split()
    df[cols] = df[cols].astype(np.uint)
    _, _, two_sided = pvalue_npy(df.calls.values, df.methylated.values, df.calls2.values, df.methylated2.values)

    t1, t2, m1, m2 = df.calls, df.calls2, df.methylated, df.methylated2

    _or =   (t1 * m2) / (t2 * m1)

    df.insert(df.shape[1], "PValue", two_sided)
    df.insert(df.shape[1], "OR", _or)

    return df

result = j.apply(fe, nb_cpu=1)
result.drop(like="methylated*|calls*")
# +--------------+-----------+-----------+------------------------+--------------------+
# | Chromosome   | Start     | End       | PValue                 | OR                 |
# | (object)     | (int32)   | (int32)   | (float64)              | (float64)          |
# |--------------+-----------+-----------+------------------------+--------------------|
# | chr21        | 5065660   | 5065810   | 0.1152954430992567     | 0.0                |
# | chr21        | 5068140   | 5068290   | 0.9999999999995065     | 0.9846955389124064 |
# | chr21        | 5128400   | 5128550   | 0.7075919513120225     | 0.4500499500499501 |
# | chr21        | 5154780   | 5154930   | 1.0                    | nan                |
# | ...          | ...       | ...       | ...                    | ...                |
# | chr21        | 46638480  | 46638630  | 0.24262295081967644    | inf                |
# | chr21        | 46639820  | 46639970  | 0.00024909594620266284 | inf                |
# | chr21        | 46692680  | 46692830  | 0.7354565511951735     | 1.1459649122807019 |
# | chr21        | 36784100  | 36784250  | 1.0                    | nan                |
# +--------------+-----------+-----------+------------------------+--------------------+
# Unstranded PyRanges object has 1,932 rows and 5 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

Thanks a lot for your reply. It seems like a decent solution but I need some more time to test it, as I'm focussing on getting a grant written these days...

Yeah, I'm writing a PR now. But with better code than here :)

I've tested that it gives the same results as your code. And it is super fast. Just wait for my PR plz

Unless you need to use Python and can't use subprocess, here's a quick CLI one-liner which sums signal from sorted BED5 files, over the genomic space where they overlap:
$ bedmap --echo --sum --delim '\t' <(bedops --merge A.bed B.bed ... N.bed) <(bedops --everything A.bed B.bed ... N.bed) > answer.bed

(Requires bash for process substitutions.)


text:
tool to visualize a collection of phylogenetic trees
I have a collection of phylogenetic trees produced from clustering DNA sequences, mostly a few nodes in each cluster/tree, with several clusters of size 1. And I would like to visualize them all together to get an idea of the relative sizes of the trees in the collection of clusters (cluster set).
E.g.
(ABC:0.0,('CL-118070':0.02010925,'CL-117615':0.00526463):0.00443474);
(DEF:0.0,('CL-117360':0.00967856,'CL-117181':0.00967856):0.00990975);
(GHI:0.0,('CL-123443':0.00123443,'CL-123443':0.00132446):0.00123435);
(JKL:0.0,('CL-134123':0.00123441,'CL-545443':0.00343446):0.00134335);
((((CL-114948:0.00908659,(CL-114989:0.00324875,CL-114959:0.00536289):0.00037038):0.00457230,(CL-114851:0.01291841,CL-114831:0.00590321):0.00077034):0.03411459,XYZ:0.00251364):0.00207594,(CL-114950:0.03074481,CL-114699:0.02256162):0.00733986,CL-116426:0.01643848);

I am thinking there must be a tool where I give it an input file with a tree per line and shows them as all in one go.
Any suggestions:
EDIT: looking at the options below:
https://en.wikipedia.org/wiki/List_of_phylogenetic_tree_visualization_software
FigTree: in this tool, one can load several trees, one per line, and there is a 'Prev/Next' menu to iterate over them, but no apparent option to display all of them at once.

Phylo.io: in this tool, one can load several trees, and display one at a time or two of them, with colouring options that compare the common nodes/branches. There isn't an apparent option to visualize many trees in one go.

PhyBin: this tool (https://cs.indiana.edu/~rrnewton/projects/phybin/) seems to cluster the sets of trees, but not visualize all of them together.
GrapeTree: this tool () can visualize only one tree at a time. It displays the last tree in the list if more than one are copy+pasted in the input text box.

There is another group of tools under Graph drawing:
https://en.wikipedia.org/wiki/Graph_drawing
Cytoscape: this tool is designed for network visualizations, which could be an option. Unsure if it visualizes newick trees, but there is are tools to convert trees to networks that could potentially work in combination with Cytoscape, e.g. https://github.com/MonsantoCo/Newick2Network

The answer to the question is to do it programmatically via Pythons ETE and use Matplotlib graphing output options to produce multiple plots per diagram. The problem is we don't know how your tree files are stored. So I can put code but I'm not sure if you could modify for your purposes

I suspect Rs ape can do the same thing but I'm not familiar with it

How would they be displayed? Are you thinking of some sort of multi-tree alignment with the trees superimposed on top of one another? Or do you just want one next to the other?

What exactly are you trying to compare? You can just count the number of sequences in each cluster for starters. With a package like ete3 you can count the number of leaves, number of nodes etc. in each tree. These kind of stats might be more informative than a visualisation depending on your interest

I recommend using densitree: https://www.cs.auckland.ac.nz/~remco/DensiTree/
And rwty: https://github.com/danlwarren/RWTY
Densitree plots the superimposed set of trees in a heatmap-like format:
And rwty, a very powerful R package that allows visualisation of treespace in an NMDS format (among other methods):

could you explain a little bit about both packages please?

One possible solution is to use ETE Toolkit Treeview, which given a set of Newick files, it will show all of them, albeit connected as shown in the screenshot:

If you think this is the answer, then we're good. I was thinking of printing unconnected trees on a single page

For comparing a large number of trees side-by-side - Dendroscope.
For combining the trees into a network - Splitstree.

Dendroscope is worth a look.  I would also look into the R phangorn package.

Dendeoscope needs a hypothesis, and interpretation is hard (for the untrained eye). Consensus phylogeny is better. Ultimately we still don't know what the OP wants to achieve

source: Stack Exchange

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

sort: votes


text:
Calculating the trimming parameter when constructing an additive phylogeny
How do you calculate the trimming parameter when constructing an additive phylogeny?
I think it has a relationship with the length of the shortest hanging edge, but I can't figure out how this works.

I would assume using the trimming parameter tool,
TreeTrimmer https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3637088/

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


text:
Profile HMM to find sequence motif using Baum-Welch instead of using a given multiple sequence alignment?
I know this can be very problem-specific -- but it seemed to me the standard approach to train an HMM to recognize a binding pattern of a protein is to use the multiple sequence alignment (MSA), instead of using Baum Welch. I suspect that this is the case because this is what I saw in the current version of HMMer3, which does not provide a baum-welch training procedure for unaligned sequence.
The thing that I saw training with Baum Welch is that once a HMM is trained, using the Viterbi can give a pattern that contains a lot of insertions in the motif, which doesn't make biological sense. This could be mitigated by using Viterbi training, however, I don't see much literature discussing this as a viable approach.
It'd be appreciated if someone could share their experience on training a profile HMM to finding motifs.

A wee caveat about the biological relevance of indels. The alignment step in HMMer3 is done with hmmalign, a local aligner that lacks a flag to change penalties as far as I known. But you can use a different aligner —but no gblocks. If done correctly, MSA indels biologically somewhat correspond to actual indel in the encoding DNA. Generally these are protein loops. And lastly, the gappy regions get really low weights the hmmbuild step, so it does not matter much. So getting too many gaps is a sign one of the sequences may be too different (non-homologous or partially in global alignment).

Could you explain why you feel B-W is a better approach? The reason we want MSAs for this is because that lets us use real data. I'm about as math challenged as a professional scientist can be, so I don't pretend to even understand what the B-W algorithm does, but as a biologist, training on real protein alignments makes intuitive sense to me. Why would you not want to go down that route?

@terdon Because HMM itself is a full probabilistic model and the typical way to train a HMM is to use BW (expectation maximization), which tries to maximize the conditional expected log likelihood. I am simply thinking about using entirely just this statistical model to tell me what is happening in the binding pattern given a set of unaligned sequences.


text:
pymol script to only select 5 chains within a distance from a reference and save the selection
I am generating 30 asymmetric units of a protein structure. Now, I want to select different combinations of 5 units which are 30 Angstroms away from the reference and save it as a different file. Here is the command I could work out till now :
fetch 3q2v
symexp sym,3q2v,(3q2v),30
select chains, bychain all within 30 of 3q2v

How to loop in and select only 5 chains at a time obtained from the setup obtained by the above code, and save it in separate files?

Firstly, I would preface that the PyMOL command line in the GUI is great, but you can run PyMOL as a regular Python library (calling the appropriate python). This makes life so much simpler for more complex manipulations. You can do really complex manipulations with it. I cannot suggest it enough.
The API version of the commands works only marginally differently, but you can do API calls in the GUI cmd.) and PyMOLWiki documentation is only slightly weaker for them (help() fills the gaps).
But if you really want to use the GUI route...
Let's get the names of the models in the selection
onames = cmd.get_object_list('chains')
print onames

Now, you want to save subsets of 5 chains (technically, objects or models). That is, you do not care about the order they are in and you don't want the same one twice (actually selections with 4 or less), so you want all combinations. Make sure to do the following
    print len(onames)
And google that number as <n> choose 5. That is how many files you will get. I am pretty sure most of these sets will be meaningful and will be useless.
In a script a loop is easy, while in the gui you have to keep it simple or use a list coprehension hackishly.
 from itertools import combinations
 for s in combinations(onames, 5): cmd.save('-'.join(s)+'.pdb', ' or '.join(s))

Obviously, you could alternatively use an enumerate method around the iterator and get an index for the file. But do note f-strings will not work.
 for i, s in enumerate(combinations(onames, 5)): cmd.save(str(i)+'.pdb', ' or '.join(s))


text:
How to analyze expression of certain group of genes in certain types of cell?
I need to analyze expression of ~ 400 genes with a certain function in the embryonic stem cells. All these genes are combined into one database with RefSeq IDs. What is a recommended workflow to quantify differential expression of these genes in RNA-Seq data of cells and visualize it?
All genes for analysis concern Human E3 Ubiquitin Ligases originally located here, hpcwebapps.cit.nih.gov/ESBL/Database/E3-ligases and the updated database is here:
https://esbl.nhlbi.nih.gov/Databases/KSBP2/Targets/Lists/E3-ligases/

By database do you mean a table? Your question is a bit broad.. have you tried anything so far?

I really apologize for it being so broad because I am a very new to bioinformatics and a little bit lost. I mean all genes I need to analyze are in a database https://hpcwebapps.cit.nih.gov/ESBL/Database/E3-ligases/

Hey no problem. Ok if you have 400+ genes, and expression values, you need to have them onto a count table, where rows are genes, and columns are your samples. If you can do some basic R, then try to run DESeq (https://www.nature.com/articles/nprot.2013.099) or edgeR (https://bioinformatics-core-shared-training.github.io/cruk-bioinf-sschool/Day3/rnaSeq_DE.pdf)

You might have a look at papers such as this one: https://www.embopress.org/doi/10.15252/msb.20188746


text:
Common QC parameters and threshold in human WGS pipeline
I am trying to summarize a list of QC-parameters and its threshold that should be applied for human WGS in general to avoid bad quality of data.
Which qc parameters for WGS are you using and why?
Thanks!

Hi and welcome to the site! Could you please [edit] your question and narrow it down a little. We don't do discussions here, we're strictly Questions and Answers, so the more specific the question, the better. What parameters are you interested in? GC content? Coverage? FASTQ quality? Insert size distribution? What will the sequences be used for (different applications may need different approaches). Please narrow the scope of your question down.

I would suggest you look at the website for the excellent MultiQC tool (https://multiqc.info/).  We use it in a clinical setting to collate all the QC statistics from our pipelines into a single report (It displays results from multiple QC programs such as FastQC, Picard, bcl2fastq, VerifyBAMid etc).  The homepage has a selection of example reports - clicking on "Whole Genome Seq" shows a report with the metrics most labs would expect from their WGS pipeline.
If you dig around in their GitHub repo (https://multiqc.info/) you will get a good idea of the default thresholds set for flagging results - there should be a config file somewhere specifying all these values.

For your negative controls, you might want to check out QC Fail (https://sequencing.qcfail.com/).  There you can see examples of how things look in QC when things have gone wrong, such as loss of quality, position-specific errors on flowcells and contamination with primer dimers.


text:
How can I export a full alignment from IGV as an image?
The IGV browser lets you export an alignment as an image (File => Save Image). However, this image only contains those reads that fit in the viewing window:

As you can see in the image above by looking at the scroll bar on the right, this is only a small subset of aligned reads. Is there any way to export an image of the aligned reads that includes all reads and not only those currently shown in the window?

Running an IGV batch script helps increase the number of reads per image, compared with File | Save image, but it does not solve the issue. To run a script, first make a script file like so with 1 command per line: snapshotDirectory ./ goto chr2:28582094-28582790 snapshot chr2_28582094_28582790.png goto chr4:8203487-8203815 snapshot chr4_8203487_8203815.png Run the script using Tools | Run batch script menu.

@TimurShtatland sorry, I have no idea what that is supposed to mean. That looks like a set of scripting commands and I have a vague recollection of having seen something like that for IGV in the past. Could you please post an answer explaining what you are suggesting?

Oops, pasted the wrong thing. Fixed this. IMO, this is not an answer, but might help partially.

@TimurShtatland yes, thanks. I'd posted my comment before your edit :). If the scripting approach does give more reads, then it's at least a partial solution and deserves an answer. Could you flesh your comment into an answer? If you'd rather not, I can post one as soon as I've tested it, but you deserve the rep.

@TimurShtatland well, this is embarrassing. Your comment made me find an old script I had which does exactly this! I had completely forgotten! Please post an answer, you should get the credit, and then I'll just edit it to add what I did to get everything (essentially, you can set maxPanelHeight to get more reads).

Thank you, done. Please feel free to edit. I don't deserve the full credit obviously, so I owe you half. :)

@TimurShtatland thanks, but that's fine. The credit belongs to Pierre whose post I'd managed to find last time I needed this and then promptly forgot about until your comment reminded me :). And I already have all the editing privileges rep can buy here, so you need it more!

Let us continue this discussion in chat.

Running an IGV batch script helps increase the number of reads per image, compared with File | Save image, but it does not solve the issue by itself.
To run a script, first make a script file like so with 1 command per line:
snapshotDirectory ./
goto chr2:28582094-28582790
maxPanelHeight 10000
snapshot chr2_28582094_28582790.png

Run the script using Tools | Run batch script menu. The importat trick is the maxPanelHeight 10000 which is what makes this include more reads than the default snapshot.
You can also convert this into a little shell script:
#!/bin/bash
mkdir -p /tmp/igv
sleep 2
echo "snapshotDirectory /tmp/igv"
sleep 2
goto chr2:28582094-28582790
sleep 2
echo "maxPanelHeight 10000"
echo "snapshot chr2_28582094_28582790.png"

Then, just open IGV, load the relevant BAM file and run the script like this to connect to IGV via telnet (credit to Pierre Lindenbaum):
sh foo.sh |  telnet 127.0.0.1 60151


text:
SRA Toolkit execution problem
I am trying to install SRA Toolkit on MacOS, following all the instructions from NCBI website. It is installed but fastq-dump cannot be executed because "cannot execute binary files". I tried "chmod +x", noah from home-brew to make binary executable on MacOS but fast-dump still cannot be executed. I attached the screenshot of my terminal window. Can anyone help?

Did you upgrade to Catalina (a bad thing due to no support for any 32bit at all)?

You can't run Linux binaries on OSX (or vice versa). Download Mac OSX binaries.


text:
Pattern mining from a genomic sequence
I need to find the following pattern from a genomic sequence
TTCGB-N15-TTCGB
B represents C, G, or T
N15 represents A, T, G or C and these four base combination should have the length of 15 bases.

To do the same, I tried grep command, but I could not specify B and N15 in the grep command line. Therefore, please help me to find the same.
For example, I have a nucleotide sequences as follows
>ATGCTTATTGCCCATTTCGTGCATGCATATGCGCATTCGCGATCGATTAGGGATAT

Expected outcome,

TTCGTGCATGCATATGCGCATTCGC

Thanks in advance.

Could you please update your question to clarify that N15 represents A, T, G or C for exactly 15 bases?

echo ">ATGCTTATTGCCCATTTCGTGCATGCATATGCGCATTCGCGATCGATTAGGGATAT" | grep -oP  'TTCG[CGT][ATGC]{1,15}TTCG[CGT]'

[CGT] looks for one occurrence of one of the three.
[ATGC]{1,15} looks for "up to 15 letter long" combinations of the four bases.

Moreover, -P makes your regex call Perl-like and is required to make this regex work. Different regex flavors behave (slightly) differently.

Dear @haci,                                                                                I have followed your command line for the test data (I had given in the example section of my question), which gives proper result. But when I tried the same with my genome file with the following command, echo $(<org15.fasta) | grep -oP  'TTCG[CGT][ATGC]{1,15}TTCG[CGT]' and I am getting output as follows,   TTCGCCGCGGCGACGTTTTCGC TTCGCATTCGG TTCGCCACCTTGTTCGG TTCGGCTTTCGC TTCGTGGCCTTCGG TTCGGCTTCTTCGT TTCGTCAAGGACCCGATTTCGG TTCGCATCTTCGC

The thing is that, the script mining TTCGB.....TTCGB, but fails to fulfil the length N15 criteria. kindly help me to fix this. Thank you.

What is your "N15 criteria" though? [ATGC]{1,15} looks for up to 15 repetitions of the four bases as you specified in your question. And as far as I can see in the output you get, this is what you are getting, 2 TTCGB stretches are connected by up to 15 bases.

N15 usually means exactly 15 of any nucleotides. So use [ACGTN]{15} . Also, -i flag for grep may need to be added if some of the sequence is lower case, for example repeat masked.

@haci N15 means any of the 4 bases, but N15 length should be 15. exactly, The output should be like this, TTCGB-ATGCTAGCATGATAG-TTCGB. The output sequence should have TTCGB at both ends, at the same time in between TTCGB 15 bases also should be there. That 15 bases can be the combination of any of the four bases, but it should have 15 bases. Sorry, if my question is not clear.

OK, then just change [ATGC]{1,15} to [ATGC]{15} as Timur Shtatland suggested. The reason for [ATGC]{1,15} was your phrase in the question: "N represents A, T, G or C up to 15 bases".

@haci I did it echo $(<org15.fasta) | grep -oP 'TTCG[CGT][ATGCN]{15}TTCG[CGT]'. It works, thank you so much @haci and @Timur Shtatland.


text:
How do researchers account for heterozygosity in genome assembly?
this is a very general question about diploid genome assembly. I am wondering how people deal with heterozygous region in a genome when assembling ?
Do they pick one of the haplotype to be in the referene and leave out the other ?
I am particularly thinking about situation similar to mammalian chromosome Y.
Do you put Y or X in the assembly ? An hybrid ?
Thanks

Generally a single allele is used. This may be arbitrary or match the major allele frequency of the species. For sex chromosomes, usually both/all are included, so for mammals both X and Y are part of the reference genome.

In order to do so, do they perform "haplotype phasing" ? So, if you have something similar to sexual chromosome on another pair (really big divergent regions), one of them will be left out, as it is not recognise as a particular pair ?

There is usually nothing special done in this regard unless it's a more finished genome (mouse, human, etc.).


text:
Find covariation of different DNA-binding protein binding between two conditions
I have ChIP-seq data of RNAPII and two other factors which we think follow RNAPII during transcription in two different conditions, and I'd like to show that the genes that lose RNAPII also lose the other two factors, and same for gains.
Do you know any methods that can calculate some sort of covariation coefficients?
I know there are methods that calculate covariation of genes between two conditions, but none for different measures on the same gene.

do you have an example dataset?

@StupidWolf Not yet, I'm still in preprocessing. I had a few issues, but I should have it soonish

Answer from @stupidwolf converted from comment:
I would try something like this, for every gene you have a logFC of RNAPII in control versus treatment. Do this for other 2 proteins, and you correlate the logFC. In very rough terms, correlation is the scaled form of covariance, so this will work.
If you have 1000 genes, do correlation between 1000 logFC from RNAPII against 1000 logFC from other proteins.
Answer from @m converted from comment:
You use a covariation stat and then perform Pearson's correlation to standardise the response between -1 to 1. An answer of >0.7 is considered robust and >0.8 highly robust. ITs an interesting question and you are right covariation is the way to answerr it


text:
How to download FASTA sequences from NCBI using the terminal?
I have following accession numbers of the 10 chromosomes of Theobroma cacao genome.
NC_030850.1
NC_030851.1
NC_030852.1
NC_030853.1
NC_030854.1
NC_030855.1
NC_030856.1
NC_030857.1
NC_030858.1
NC_030859.1

I need to download these FASTA files using the terminal because I'm working on a server. I'll be grateful if someone can help me.
Thank you.

Alternatively, you can use the NCBI Entrez Direct UNIX E-utilities
Basically, you have to download the install file here: https://www.ncbi.nlm.nih.gov/books/NBK179288/bin/install-edirect.sh
In the terminal, install it using: source ./install-edirect.sh
Then, you can download your sequence by doing:
esearch -db nucleotide -query "NC_030850.1" | efetch -format fasta > NC_030850.1.fasta
And you should find your fasta sequence downloaded.
As you have several sequences to download, I think it will be quite easy to add this command into a little bash script to process all of them.
More readings about that: http://bioinformatics.cvr.ac.uk/blog/ncbi-entrez-direct-unix-e-utilities/ and https://www.ncbi.nlm.nih.gov/books/NBK179288/

Thank you very much for all the support.

@MudithMMBc, you're welcome ;)

@MudithMMBc Try to get a NCBI_API_KEY for batch downloading with GNU parallel at a max rate of 10 sequences/s. Then save accessions in one file and run parallel -a accession.txt -j8 'efetch -db nuccore -id {} -format fasta >{}.fa'. Adjust -j8 to -j3 if you do not have a key.

Just wget or curl each as https://www.ncbi.nlm.nih.gov/nuccore/NC_030853.1?report=fasta&log$=seqview&format=text...

Thank you very much.

This doesn’t seem to work (any longer?). For instance, curl -L 'https://www.ncbi.nlm.nih.gov/nuccore/NC_030853.1?report=fasta&log$=seqview&format=text' downloads a HTML page with an (invisible) form on it, not a FASTA file.

Downloading entire genomes
The best way to download FASTA sequences for an entire genome is to search for the genome, for example Theobroma cacao genome in the NCBI Assembly portal and use the big blue Download button.

Downloading individual chromosomes
For a given assembly, if you want to download the FASTA sequences for a bunch of chromosomes, you can do that by going to the Genomes FTP path highlighted in the screenshot:

Once you are in the Genomes FTP path, you can navigate to the FASTA folder as follows:
Theobroma FTP directory >> Assembly Structure >> Primary Assembly >> Assembled Chromosomes >> FASTA.
Downloading a few sequences
For this, you can use Entrez Direct as mentioned by @dc37. Why not always use Entrez Direct? While it is fine for a small number of sequences, it can be slow to download a large number of sequences. Entrez Direct by default will download uncompressed data so you will end up spending more time downloading a larger file instead of downloading a smaller, compressed file from FTP more quickly.
If you were to use Entrez Direct for this purpose, I'd not bother with a bash script and use epost to first post the entire list of accessions and then pipe it to efetch as shown below:
epost -db nuccore -input <file_with_accs> -format acc \
  | efetch -format fasta > sequences.fasta

Thank you very much for all the support.

Go to the main FASTA sequence page, view source, and get the "ncbi_uidlist" value.  In this case, it's 1061361601
Use the ID in the query string to download directly:
https://www.ncbi.nlm.nih.gov/sviewer/viewer.fcgi?id=1061361601&db=nuccore&report=fasta&retmode=text&withmarkup=on&tool=portal&log$=seqview&maxdownloadsize=100000000

You can also get this link directly, by using Chrome's developer tools (F12), viewing the Network tab, then loading the page in (1).

(Sorry, I can't post a comment with my reputation score).
@vkkodali, interesting command I did not know about it ;). However, your command is downloading all sequences from the input file into a single fasta file. I believe that if you have a big bunch of sequences, it could be a little bit tricky after to manipulate that kind of file. Just by curiosity, do you have a trick to create one sequence per file using efetch ?

If you have a big bunch of sequences, you will end up with a lot of files. Whether you want a large number of files or just one file is, I guess, a personal choice. A multifasta file is fairly standard though. I don't think you can create individual files for each sequence using epost and efetch; you will have to either use a bash script or postprocess the efetch output using the unix tool split.

Ok, I understand it. Thanks for your thoughts ;)


text:
How to interpret Mendelian randomization results?
I did Mendelian randomization using this software:
https://cran.r-project.org/web/packages/MendelianRandomization/vignettes/Vignette_MR.pdf
library(MendelianRandomization)
f=read.table("246LDout272Biobank_Retina.txt", header=T)
> head(f)
          rs exposure.beta exposure.se outcome.beta outcome.se
1  rs1029830      0.723525  0.03026430  0.066715400  0.0359278
2  rs1029832      0.723785  0.03029603  0.064105600  0.0359021
3 rs11078374     -0.411789  0.04189295 -0.000376929  0.0406439
4 rs11078382      0.882549  0.14275799 -0.197074000  0.1247720
5  rs1124961     -0.333763  0.05589377 -0.075468600  0.0576012
6  rs1135237     -0.316831  0.05552530 -0.074086200  0.0573111

MRInputObject <- mr_input(bx = f$exposure.beta,bxse = f$exposure.se,by = f$outcome.beta,byse = f$outcome.se)
EggerObject <- mr_egger(MRInputObject,robust = FALSE,penalized = FALSE,correl = FALSE,distribution = "normal",alpha = 0.05)

MR-Egger method
(variants uncorrelated, random-effect model)

Number of Variants =  246

------------------------------------------------------------------
  Method     Estimate Std Error  95% CI       p-value
  MR-Egger    0.115     0.016  0.084, 0.146   3.28e-13
 (intercept) -0.010     0.008 -0.025, 0.006   0.226
------------------------------------------------------------------
Residual Standard Error :  1.106
Heterogeneity test statistic = 298.4508 on 244 degrees of freedom,     (p-value = 0.0099)
 I^2_GX statistic: 97.2%

Let's say that my exposure is called Retina and my outcome is called Biobank. Can someone please help me interpret these results in terms of 'horizontal pleiotropy' as mentioned here: https://www.ncbi.nlm.nih.gov/pubmed/29771313

It looks like you have linkage as opposed to linkage disequilibrium, therefore horizontal pleiotropy is an explanation of the calculation.
The calculation appears to be a GLMM around the beta distribution, the residual looks low (which is good). The heterogeneity statistic is significant, you would need to assess whether this is an index of heterozygosity, but it will relate to expectation of Mendellian inhertiance based around the Hardy Weinberg equilibrium.


text:
working with cutadapt
I'm working with ion torrent data where I apply the program cutadapt.
I'm analyzing ITS seq data, as well as Matk.
When I'm using cutadapt, I search the information for this genes, and built new files accordingly to his information. For instance, when I search for ITS, I generate a file about this gene, and the same applies to the Matk gene.
However, more recently I randomly picked a fasta sequence from ITS, blast it and received the information about the matk gene. Is this supposed to occur?
Thank you for the help and so sorry for the naive question.

How do you know that you didn't just mix up your Matk and ITS data? Not sure that I understand how you "generate a file".

Answer from @m, converted from comment:
It's definitely not supposed to happen. There appears to be a deep rooted classification bug, where matk is misclassified as ITS. Unpicking this will be hard.


text:
Method to detect genome doubling
CNV is usually detected using either read depth or SNP VAF. However, neither of these methods can detect CNV in case of genome doubling (assuming genome doubling is a kind of genome scale CNV comparing to normal tissue which is diplod). We thought of other platform such as SNP array and we did not think they can detect genome doubling as well. Can any one share some comments on the method for such event? Thanks in advance.

Are you talking about SNVs or CNVs? I think you are asking about CNV detection in species that have undergone a while genome duplication event, for example teleost fish. Is that what you mean? Or are you working on a novel genome assembly and want to see if the genome of this species has been duplicated? Please [edit] your question and clarify. If you give us more detail, we might be able to answer better.

In a tissue or cell line that used to be diploid and has undergone genome wide doubling, you expect more variants having observed allele frequencies that would also occur in tetraploid organism. That is, observed VAF distribution will have (possibly small) peaks at 0.25 and 0.75, in addition to the diploid peaks at 0, 0.5 and 1. Compare the distributions in the tissue or cell line of interest to those of the control, that is, known diploid. Also if it is a doubled male genome, the sex chromosomes will show a new peak at AF 0.5, that is, heterozygous X or Y variants. This is just an idea, I have never tried to look for such events.

The following may work, assuming you have an Illumina genome sequence library with decent coverage, preferably amplification-free:
Use kat to produce a kmer spectrum plot of the reads (kmer abundance vs. kmer frequency). When you look at, for example, the 21-mer spectrum of the reads, you would expect to see a peak in the spectrum where the kmer frequency corresponds with the estimated coverage of the genome.  For instance, if your reads would be expected to provide 10x coverage of the original diploid genome, then you would expect a peak in the kmer spectrum at a kmer frequency around 10.
If the genome had been doubled, and each of the two copies had undergone a number of mutations, then you would also expect a second peak to be found with half the coverage of the first (at a kmer frequency of 5 in the example).  The relationship between the heights of the two peaks would reflect the degree to which diploidisation had occurred post-doubling.
Also, I seem to remember something like a pachytene chromosome spread being used with an optical microscope.


text:
Replacing some strings in R or terminal
I have a .txt file; I have 139 unique sample IDs in sampleID column
> head(mutations)
                                                                              sampleID chr     pos ref mut
4 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1  920552   G   A
5 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 1350208   G   A
6 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 1447367   T   G
7 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 1909310   G   A
8 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 2172675   G   C
9 LP2000107-DNA_A01_vs_LP2000102-DNA_A01__pv.1.6__rg.grch37_g1k__al.bwa_mem__.snp.pass   1 2846623   C   T
>

> dim(mutations)
[1] 3574407       5
>

I want to change name of sample from 1 to 139 in R or Linux terminal but I don't know how to do that, any help please?
Thanks

I believe that
uniqsampleID = 1:length(unique(mutations$sampleID))
names(uniqsampleID) = unique(mutations$sampleID)
mutations$sampleID = uniqsampleID[mutations$sampleID]

should work, assuming mutations is a data frame.

Explanation:
1:N creates a vector of length N, containing 1, 2, ..., N. So using 1:length(X) will do this for each item in X, so 1:length(unique(mutations$sampleID)) should give you a vector of 1:139.
Then names(uniqsampleID) = unique(mutations$sampleID) assigns to each of these numbers one of these sample IDs, which allows you to retrieve the sample number with uniqsampleID[mutations$sampleID].

Sorry that gives 3574407 sample IDs while I only have 139 samples for each I have 25 chromosome and a lot of mutations; I mean sampleID must contains maximum 139 unique numbers

Ah I see! I'll change it.

Yessssss worked thanks; Sorry in chr column I have the number of chromosome 1:23 plus to X and Y; How I can add chr prefix to this column? I mean instead of 1 I should have chr1

Technically should be another question, but paste0() can handle that. It would probably be mutations$chr = paste0("chr", mutations$chr)


text:
How I loop over this code
I have 139 samples; In R for each sample I run below codes
for example for sample 1
sample_1 = whichSignatures(tumor.ref = sigs.input,
                           signatures.ref = signatures.nature2013,
                           sample.id = 1,
                           contexts.needed = TRUE,
                           tri.counts.method  = 'genome')

a1=(sample_1[[1]])

> class(sample_1)
[1] "list"

> a1=(sample_1[[1]])

> head(a1)
  Signature.1A Signature.1B Signature.2 Signature.3 Signature.4 Signature.5 Signature.6 Signature.7 Signature.8
1    0.5370653            0           0           0           0           0           0           0           0
  Signature.9 Signature.10 Signature.11 Signature.12 Signature.13 Signature.14 Signature.15 Signature.16 Signature.17
1           0            0            0            0            0            0            0            0    0.3383947
  Signature.18 Signature.19 Signature.20 Signature.21 Signature.R1 Signature.R2 Signature.R3 Signature.U1 Signature.U2
1    0.1123409            0            0            0            0            0            0            0            0
> dim(a1)
[1]  1 27

and for sample 2
sample_2 = whichSignatures(tumor.ref = sigs.input,
                           signatures.ref = signatures.nature2013,
                           sample.id = 2,
                           contexts.needed = TRUE,
                           tri.counts.method  = 'genome')

a2=(sample_2[[1]])

And finally binding a1,a2,...a139 by
sig=rbind(a1,a2)

> print(sig)
  Signature.1A Signature.1B Signature.2 Signature.3 Signature.4 Signature.5 Signature.6 Signature.7 Signature.8
1    0.5370653            0           0           0           0           0           0           0           0
2    0.7162292            0           0           0           0           0           0           0           0
  Signature.9 Signature.10 Signature.11 Signature.12 Signature.13 Signature.14 Signature.15 Signature.16 Signature.17
1           0            0            0            0            0            0            0            0    0.3383947
2           0            0            0            0            0            0            0            0    0.2837708
  Signature.18 Signature.19 Signature.20 Signature.21 Signature.R1 Signature.R2 Signature.R3 Signature.U1 Signature.U2
1    0.1123409            0            0            0            0            0            0            0            0
2    0.0000000            0            0            0            0            0            0            0            0
>

sigs.input is a common R object here
Doing these manually is too painful; Can you please help me how I do that automatically in R? I am not sure but something like looping or whatever you guess possible
Thanks

list.of.results <- lapply(seq(INSERT_TOTAL_NUMER_OF_SAMPLES_HERE),function(sample.id){
whichSignatures(tumor.ref = sigs.input,
                           signatures.ref = signatures.nature2013,
                           sample.id = sample.id,
                           contexts.needed = TRUE,
                           tri.counts.method  = 'genome')[[1]]
})

sig <- do.call(rbind,list.of.results)

do you know how to do this in dplyr ?


text:
Simulating phenotype with the 1000 Genomes Project
I'm looking for a way to simulate phenotypes against a real SNP data source, such as the 1000 Genomes. It must be free for commercial purpose (Eg.: MIT license). Any recommendation?
I'm trying to use the GCTA64, but I couldn't get it working. The documentation doesn't help much.
At the end of the day, I want:
1 - Simulate Case/Control and/or quantitative phenotype
2 - Link it with a real SNP dataset (eg: 1000 Genomes)
3 - Conduct GWAS analysis using Plink and/or Hail.
Looking Hail's and Plink's tutorials for GWAS, I realised both use simulated phenotypes from real SNP data sources (1000 Genomes and HapMap, respectively), but how they created the datasets are beyond of the scope of the tutorials, thus not reported.
As mentioned before, I've tried gcta64, but no success. Here's what I've tried:
1 - Downloaded 1000 Genome sample from Plink page:
Entire dataset as a single .tar.gz (1.12 GB) (A2 allele major, not ref, on chr3 before 15 Oct 2017)
2 - Tried to generate the simulate data by:
./gcta64  --bfile 1kGenomesP1/1kg_phase1_all/1kg_phase1_all --simu-qt  --simu-causal-loci causal.snplist  --simu-hsq 0.5 --simu-rep 3  --keep test
.indi.list --out 1kg_phase1_all

Error: Error: --keep test.indi.list not found.
Question: What should be the files:causal.snplist, test.indi.list, teest? I know the property name suggests what should be, but there's no example of file structure nowhere... Any practical example or tutorial?
Btw - I apologise in advance if this is a too trivial question. I'm quite new at it. Appreciate your patience and help :)

I managed to get it done without the gcta64.
I have written a Jupyter Notebook with the steps: https://github.com/bambrozio/bioinformatics/blob/master/utils/1k-genomes-phenotype-simulator.ipynb

You might try ldscsim which was developed by Nik Baya for this purpose.
https://hail.is/docs/0.2/experimental/ldscsim.html
import hail as hl

mt = hl.import_vcf('path/to/your/data.vcf')
mt.write('data.mt') # store in Hail's native format for speed
mt = hl.read_matrix_table('data.mt')
mt = hl.experimental.ldscsim(
    mt, mt.GT, h2=0.5
)

ldscsim can generate phenotypes under a number of different models.
For the GWAS tutorial, AFAIK, we just simulated betas for each site, then computed the phenotype from that. If you ask on https://discuss.hail.is and tag Tim Poterba he can probably give you the precise answer.
def logistic(x):
    return 1 / (1 + hl.exp(-x))

mt = mt.annotate_rows(beta = hl.rand_norm(0, 1))
mt = mt.annotate_cols(
    is_case = 0.5 < logistic(hl.agg.sum(beta * mt.GT.n_alt_alleles()))
)


text:
Error with BWA Mem input having multiple fastq files using cat and process substitution
Running BWA Mem on a number of paired end fastq files using process substitution on the inputs results in this error:
[mem_sam_pe] paired reads have different names: "HWI-D00360:5:H814YADXX:1:2111:20284:33068", "HWI-D00360:5:H814YADXX:1:1213:17580:34847"

The command:
bwa mem -M ref.fa "<(cat U0a_CGATGT_L001_R1_004.fastq U0a_CGATGT_L001_R1_005.fastq U0a_CGATGT_L002_R1_001.fastq U0a_CGATGT_L002_R1_002.fastq U0a_CGATGT_L002_R1_003.fastq U0a_CGATGT_L001_R1_001.fastq U0a_CGATGT_L001_R1_002.fastq U0a_CGATGT_L002_R1_004.fastq U0a_CGATGT_L001_R1_003.fastq U0a_CGATGT_L002_R1_005.fastq)" "<(cat U0a_CGATGT_L001_R2_003.fastq U0a_CGATGT_L002_R2_005.fastq U0a_CGATGT_L001_R2_004.fastq U0a_CGATGT_L001_R2_005.fastq U0a_CGATGT_L002_R2_001.fastq U0a_CGATGT_L002_R2_002.fastq U0a_CGATGT_L001_R2_001.fastq U0a_CGATGT_L002_R2_003.fastq U0a_CGATGT_L001_R2_002.fastq U0a_CGATGT_L002_R2_004.fastq)" > out.sam

However the read (there were several) that threw this error does seem to be properly paired:
grep 'HWI-D00360:5:H814YADXX:1:2111:20284:33068' *.fastq
U0a_CGATGT_L001_R1_004.fastq:@HWI-D00360:5:H814YADXX:1:2111:20284:33068 1:N:0:CGATGT
U0a_CGATGT_L001_R2_004.fastq:@HWI-D00360:5:H814YADXX:1:2111:20284:33068 2:N:0:CGATGT

I can see that the concatenation of the input files is happening in a different order, so this is the next thing I'm going to try to fix, but has anyone else run into this?
(I probably just answered my own question but I'll post it anyway, that may or may not be the issue)

In 2019 you almost certainly don't want to be using the -M flag. See bwa #26, #96, #188, et al.

Ok thanks @JohnMarshall I'll give it a try without. I probably had it there for the older Picard or GATK but I've since updated both

You are correct that the concatenation of input read files will confuse BWA. BWA depends on the input reads being in the same order in R1 and R2 files for paired-end alignment.
From bwa docs: If mates.fq is present, this command assumes the i-th read in reads.fq and the i-th read in mates.fq constitute a read pair. BWA will not search through FASTQ files trying to find matching mate IDs.
You will notice that your initial error message is reporting two different read IDs that are improperly paired. Presumably the mispaired read corresponds to the first read in the out-of-order file.

AS both the OP and @MaximilianPress suggested, the paired-end reads should be provided as 2 separate arguments to bwa. This is true whether the reads are in files or streamed. For the sake of completeness, I provide the corrected command here.
bwa mem -M ref_file.fa <(cat A.read1.fq B.read1.fq ...) <(cat A.read2.fq B.read2.fq ...)
The order of the files is important: if the order is A, B for read1, then it should be A, B and not B, A for read2.
EDIT:
(Thanks to @JohnMarshall for pointing this out.)
Avoid double quotes around process substitution. Quote the file names (if needed). The quotes around process substitution, as apparently intended, work only on some systems/versions1, and seem like an undocumented feature to me.
DON'T:
bwa mem -M ref_file.fa "<(cat A.read1.fq B.read1.fq ...)" "<(cat A.read2.fq B.read2.fq ...)"

1
Oddly enough, process substitution inside double quotes works at least in the case I tested before posting: bwa mem -M test1.fa "<(cat test_assembler_tiny_1.1.fastq)" "<(cat test_assembler_tiny_1.2.fastq)" This contradicts the answer here: https://unix.stackexchange.com/q/486019 . I used bwa 0.7.17-r1188, with both GNU bash, version 3.2.57(1)-release (x86_64-apple-darwin17) and zsh 5.7.1 (x86_64-apple-darwin17.7.0), on macOS High Sierra 10.13.6. From OP's error msg, it seems it worked for @glyph too.

Quotes around ordinary filenames are optional, but here you've put quotes around the shell redirection metacharacters. This leads to <(cat …): No such file or directory. Don't quote the <() characters if you want them to be interpreted by the shell.


text:
How do I increase the minimum dot size in Seurat's DotPlot function?
I want to use the DotPlot function from Seurat v3 to visualise the expression of some genes across clusters. However when the expression of a gene is zero or very low, the dot size is so small that it is not clearly visible when printed on paper. Please is there a possibility to increase the minimum dot size in the DotPlot function to make the dot sizes more visible when printed?
For example, I would like to have a minimum dot size set to be like

from this publication
Rather than:

from this other publication
Thank you in advance for your helpful hint.

You can set the scale using the dot.scale option. This will affect all the points, not only the 0s.

@fra. Thank you but when I increase the dot.scale parameter,only the bigger points really change. The smaller points change only when the dot.scale value is really high and the rest of the image now looks unappealing.

It would be much easier to answer your question if you provided a reproducible example

Had also been asked on https://github.com/satijalab/seurat/issues/2297 .

Did you try to use DotPlot(..., scale.by = "size")?
In contrast to the default scale.by= "radius", this will link the area (==2*pi*r^2), not the radius, of the circles to the fraction of cells expressing the feature.
This corresponds much better to our perception of size and will make differences in low values easier to see.
Reading ?Seurat::DotPlot the scale.min parameter looked promising but looking at the code it seems to censor the data as well.
Since Seurat's plotting functionality is based on ggplot2 you can also adjust the color scale by simply adding scale_fill_viridis() etc. to the returned plot. This might also work for size. Try something like:
DotPlot(...) + scale_size(range = c(5, 10)) # will likely warn about supplying the same scale twice

Apart from this, Seurat's plotting system is not very hackable and I find it much easier to extract the relevant data and plot them myself with ggplot2.
Yet another comment: Your plot with the strong differences looks much more convincing to me wrt. to the marker property of these genes than the cited plot.

You can define your breaks with:
p + scale_size(breaks = c(0, 25, 50, 75))

You can adjust the size of the points with:
p + scale_size(range = c(1,8))     # Note: c(1,6) is the default

For more info, see ??scale_size


text:
There's a faster way to download 1000 Genomes phase 3?
I'm trying to download the 1000 Genomes (phase 3) through FTP, but it's taking ages.
I saw the documentation about using Aspera (fasp protocol), but the instructions don't work: ERR [ascp] SSH authentication failed.
Eg:
ascp -i /home/ibmuser/.aspera/connect/etc/asperaweb_id_dsa.putty -Tr -Q -l 100M \
     -P33001 -L- fasp-g1k@fasp.1000genomes.ebi.ac.uk:vol1/ftp/release/20130502/ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz.tbi ./

Via Aspera Desktop client I got:
SSH_MSG_DISCONNECT: 2 Too many authentication failures
Thanks!

You are using the command for older versions. In the page you link to, they give two commands, one for the newer and one for the older versions of aspera. You are using the one for the older. The only difference is the keyfile, but what you want is:
ascp -i /home/ibmuser/.aspera/connect/etc/asperaweb_id_dsa.openssh -Tr -Q -l 100M \
 -P33001 -L- fasp-g1k@fasp.1000genomes.ebi.ac.uk:vol1/ftp/release/20130502/ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz.tbi ./

At least, I just tried this a second ago with the aspera bundle I downloaded today and that's the command that worked: asperaweb_id_dsa.openssh and not asperaweb_id_dsa.putty.

Yep! I was getting back here to post the answer as I realised the same lol. Accepted your answer in this case. Thanks! :)

Btw - I've created a script to download everything: https://github.com/bambrozio/bioinformatics/tree/master/utils

@BrunoAmbrozio nice! But note that you forgot the sex chromosomes. And also the various alt ones if any are present.


text:
GenomicRanges: get nearest neighbor distances for random genes using a for loop
I want to compute the nearest neighbour distances using a for loop. I can do this with a random subset by doing this:
df <-  read.table("tmp.txt", header=F)
colnames(df) <- c("gene", "seqnames", "start", "stop", "LOCID")
gr <- makeGRangesFromDataFrame(df, keep.extra.columns=TRUE)

# genes of interest in a list
gene.list <- read.table("careplastic_diffs.txt")

# filters the GenomicRanges object with genes of interest
sub.gr <- gr[(elementMetadata(gr)[, "gene"] %in% gene.list$V1)]

# compute distances
dists <- distanceToNearest(x = sub.gr)

mcols(dists)[,1] # contains a list of distances

I now want to essentially do this 5 times, using the same logic as above:
## create a rand dataframe, each column is a random list of 50 genes
rand.df <- data.frame(replicate(5, sample(gene.list$V1, 50)))

# loop through columns and gets distances and prints to a file
for(i in 1:ncol(rand.df)){
  rand.gr[i] <- gr[(elementMetadata(gr)[,"gene"] %in% rand.df[,i])]
  dists[i] <- distanceToNearest(x = rand.gr[i])
  output <- print(mcols(dists[i][,1]))
}

But I can't seem to make this work.

Error: subscript contains out-of-bounds indices In addition: Warning
  message: In recycleSingleBracketReplacementValue(value, x, nsbs) :
  number of values supplied is not a sub-multiple of the number of
  values to be replaced

It's probably something really silly but any help would be appreciated ...

Is it because rand.gr has not been initialised, so when you ask to set rand.gr[i] it can't because there isn't one yet?
If you set i to a value (e.g. 1), and run the lines inside your loop, which of the lines produces the error?
If you don't need to retain all the random genomic ranges at the end, you could just leave out some of the [i]s, something like this:
for(i in 1:ncol(rand.df)){
  rand.gr <- gr[(elementMetadata(gr)[,"gene"] %in% rand.df[,i])]
  dists <- distanceToNearest(x = rand.gr)
  output <- print(mcols(dists[,1]))
}

Isn't output being overwritten with each iteration?

Thanks, I do need a file with essentially columns of distance data (each column corresponding to a column in rand.df)

this seemed to create another problem: https://www.dropbox.com/s/jl03i9nm87vwbxz/Screenshot%202019-11-08%20at%2013.21.04.png?dl=0

haci is correct, ouput is overwritten on each time round the loop.

You could do this much more easily and much faster on the command line, starting with an unsorted BED file of genes called genes.unsorted.bed:
$ sort-bed genes.unsorted.bed | shuf - | head -5 - | sort-bed - | closest-features --closest --dist --no-overlaps - > answer.bed

Replace the 5 in head -5 with however many random genes you want. Remove --no-overlaps if you want to allow consideration of overlapping genes (which returns a distance of zero).
Output columns are delimited with |. This can be replaced with another delimiter, but the output can be brought into R easily with read.table or similar.
Ref: https://bedops.readthedocs.io/en/latest/content/reference/set-operations/closest-features.html

i thought closest-features requires a reference and an input file?


text:
HVG and variation coefficient
I'm looking for highly variable genes from a data set. I have calculated the coefficient of variation of the expressed genes of my dataset and inserted them in a histogram graph. Now I should select only the highly variable genes. Should I select genes with a coefficient of variation >1, or how?
thank you for your help,
regards

Heuristically, I agree with @Jonathan Moore's answer. However, you could try to take a more disciplined approach and model the distribution as a mixture of two distributions, and infer the parameters of the two distributions from that, which will allow you to make somewhat more informed decisions. A basic intro to mixture modeling in R is here.
However, a worry of mine when I see data such as this is that you have some confounder that leads to that spike at CV ~ 1. I would actually plot the mean against the SD of the genes in question, and make sure that it is not due to observations of e.g. one transcript across samples when in most samples there are zero transcripts (see last section here). This would almost certainly be noise, and then what you would be enriching for is genes with very low expression, leading to a very small denominator for the CV calculation. If you are pseudocounting the expression (for example adding 1 to all values to avoid dividing by zero), I would try a few different values for the pseudocount to ensure that the value of the pseudocount is not significantly affecting your CV distribution that you show here. That would be a major warning sign of artifacts.
Therefore I would also suggest setting some sort of cutoff for minimal expression of genes in that distribution. Maybe an average of 1 transcript per observation (sample [bulk RNA-seq] or cell [scRNA-seq]) might be a reasonable cutoff for genes to consider.
I would do this data cleaning and filtering before making any decisions about this CV distribution.

thank you very much for the suggestion. i have already made a cutoff, but probably it was not enough

There seems to be a dividing line around CoV=0.9.  Perhaps this is a good place to start for a cutoff.


text:
how to calculate Pairwise Alignment Scores for blosum62
I have a bioinformatics exam coming up. I can understand the difference between the global algorithm and the local algorithm, but I have a problem with gap opening penalty and gap extension penalty. What is the difference between gap opening penalty and gap extension?
Can anybody help me in working out how to use the gap opening penalty and gap extension penalty for blosum62?
There are a lot of questions like the one below:
Compute the score of the following sequence alignment given the BLOSUM62 matrix below and gap opening penalty (GOP) = 12, and gap extension penalty (GEP) = 2.
SDRVIKAAIFDPIQPDF---G-----PVYFGLGHVH
RDLVERLFILDMI-PGLIKAGDSFPIPVALMINHIF

   A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V
A  4 -1 -2 -2  0 -1 -1  0 -2 -1 -1 -1 -1 -2 -1  1  0 -3 -2  0
R -1  5  0 -2 -3  1  0 -2  0 -3 -2  2 -1 -3 -2 -1 -1 -3 -2 -3
N -2  0  6  1 -3  0  0  0  1 -3 -3  0 -2 -3 -2  1  0 -4 -2 -3
D -2 -2  1  6 -3  0  2 -1 -1 -3 -4 -1 -3 -3 -1  0 -1 -4 -3 -3
C  0 -3 -3 -3  9 -3 -4 -3 -3 -1 -1 -3 -1 -2 -3 -1 -1 -2 -2 -1
Q -1  1  0  0 -3  5  2 -2  0 -3 -2  1  0 -3 -1  0 -1 -2 -1 -2
E -1  0  0  2 -4  2  5 -2  0 -3 -3  1 -2 -3 -1  0 -1 -3 -2 -2
G  0 -2  0 -1 -3 -2 -2  6 -2 -4 -4 -2 -3 -3 -2  0 -2 -2 -3 -3
H -2  0  1 -1 -3  0  0 -2  8 -3 -3 -1 -2 -1 -2 -1 -2 -2  2 -3
I -1 -3 -3 -3 -1 -3 -3 -4 -3  4  2 -3  1  0 -3 -2 -1 -3 -1  3
L -1 -2 -3 -4 -1 -2 -3 -4 -3  2  4 -2  2  0 -3 -2 -1 -2 -1  1
K -1  2  0 -1 -3  1  1 -2 -1 -3 -2  5 -1 -3 -1  0 -1 -3 -2 -2
M -1 -1 -2 -3 -1  0 -2 -3 -2  1  2 -1  5  0 -2 -1 -1 -1 -1  1
F -2 -3 -3 -3 -2 -3 -3 -3 -1  0  0 -3  0  6 -4 -2 -2  1  3 -1
P -1 -2 -2 -1 -3 -1 -1 -2 -2 -3 -3 -1 -2 -4  7 -1 -1 -4 -3 -2
S  1 -1  1  0 -1  0  0  0 -1 -2 -2  0 -1 -2 -1  4  1 -3 -2 -2
T  0 -1  0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1  1  5 -2 -2  0
W -3 -3 -4 -4 -2 -2 -3 -2 -2 -3 -2 -3 -1  1 -4 -3 -2 11  2 -3
Y -2 -2 -2 -3 -2 -1 -2 -3  2 -1 -1 -2 -1  3 -3 -2 -2  2  7 -1
V  0 -3 -3 -3 -1 -2 -2 -3 -3  3  1 -2  1 -1 -2 -2  0 -3 -1  4

so could you help me please about solution?

The lower matrix is the blosom62. The two sequences are the pairwise alignment... easy. The gap penalties are described by here... https://en.m.wikipedia.org/wiki/Gap_penalty ... I don't have code for this and I'm not a fan of optical analysis, so I'm not going to add up the numbers

my function align not using matrix... i have implemented one algo that can fail...possible it is better one good eye too

Your alignment is

SDRVIKAAIFDPIQPDF---G-----PVYFGLGHVH
RDLVERLFILDMI-PGLIKAGDSFPIPVALMINHIF

For each position in the alignment you calculate the score for that alignment. For position 1 we'd look up S vs R in the matrix and find a score of -1. You continue doing this until you hit the first -, which is not in the matrix. The first - is a gapopening, each consequent - in a series of -'s counts as a gap extension, instead of an opening.
---G--
IKAGDS

gives
-12 -2 -2 6 -12 -2

The problem with the example is a gap penalty of 12 and 2 means the position of G is critical to the output, regardless of its location as a biological homologue. If G is at the beginning or end of the indel the score drops considerably, ie only a single gap opening penalty.

That is a matter of alignment score maximisation, the author's question regards calculating the score of the given alignment with the given BLOSUM matrix, and how to incorporate the gap penalties into this.

SDRVIKAAIFDPIQPDF---G-----PVY------FGLGHVH
RDLVERLFILDMI-PGLIKAGDSFPIPVALMINHIF------ has one uguagliance more


text:
What is a samtools mpileup reference skip?
The samtools documentation for mpileup states:

At this column, a dot stands for a match to the reference base on the forward strand, a comma for a match on the reverse strand, a '>' or '<' for a reference skip
  ...
  Similarly, a pattern `-[0-9]+[ACGTNacgtn]+' represents a deletion from the reference.

What is a reference skip?

That documentation has been expanded recently (see PR #1055) and now describes the characters seen in this column in more detail:
Forward   Reverse                    Meaning
---------------------------------------------------------------
 . dot    , comma   Base matches the reference base
 ACGTN     acgtn    Base is a mismatch to the reference base
   >         <      Reference skip (due to CIGAR "N")
   *         *      Deletion of the reference base (CIGAR "D")

(And so on — see the pull request or doc/samtools-mpileup.1 on the develop branch for the complete pileup format description. The documentation on that web page will be updated when this is released in a new samtools release.)
As noted in https://bioinformatics.stackexchange.com/a/421/134, SAM has always had separate D (deletion from the reference) and N (skipped region from the reference) operators. These are treated similarly computationally but emitted intentionally by aligners to distinguish introns and the like from other deletions.

Skipped references are similar to deletions, but the different symbols used are indicating that it's an expected deletion. The most common use for a skipped reference is when excluding intronic sequence from cDNA reads that are being matched to a genome reference.
Note that the skipped reference bases are added by the mapping program, which generates the SAM/BAM/CRAM format and decides on whether to call a deletion or a skipped base.

How would samtools recognize that? Surely it doesn't mode splice sites! Does this mean we only see a reference skip if using a particular samtools flag to tell it we're feeding it cDNA? Otherwise, won't it just see a deletion?

Samtools is a set of tools to operate on SAM/BAM/CRAM files; it doesn't generate them in the first place. I've added clarification to my answer.

thanks, I know what they are. What I'm wondering is how samtools will differentiate between a deletion and a skipped reference. Both will appear as a gapped alignment, presumably, so what's the difference between one and the other in the bam/sam file?

The person who asked the question already mentioned that:
"a '>' or '<' for a reference skip... a pattern `-[0-9]+[ACGTNacgtn]+' represents a deletion from the reference"

I know. I guess I must be missing something really obvious, but what I'm asking is how samtools differentiates between a normal deletion and a reference skip. Since a reference skip is just a deletion, I don't see how samtools would know to mark it specially. Presumably, the necessary information is somehow encoded in the bam, but what is that? How can samtools know whether this deletion is an expected one, as you put it?

Repeating myself a bit here... samtools doesn't know this. That information is provided (or not) by the mapper.

So it is somehow encoded in the bam file, thanks. That's what I was wondering about.

@terdon: The SAM specification is at https://samtools.github.io/hts-specs/SAMv1.pdf. CIGAR operators are described in §1.4 or you can search for “intron”.


text:
Gene Ontology Specific Fold Changes
I am currently doing a project involving gene rank correlations.  How we rank our genes is based on the absolute values of the log2 fold change.  I.e. a gene which is doubled with respect to the control condition is listed as log2FC=1 a gene at 4xbaseline would be log2FC=2 and a gene at half baseline would be log2FcC=-1.  Since we are more interested simply in if there is an effect i.e. we are hypothesis free about whether or not gene under expression or over expression should be considered important we took the absolute values of these genes so hypothetically a gene at half the expression level of baseline are ranked the same as genes at double the expression level of baseline.
My concern is that this list is across all gene types, i.e. we are lumping kinases in with ion channels and axonal binding etc etc.  Would certain gene ontologies have different log fold changes simply due to a smaller range of raw expression?  And therefore should be normalized within their ontology?

Ouch, that's going to be tricky. Especially considering that the majority of proteins have multiple, and often very separate, molecular function annotations. You won't find it easy to separate them into clean bins based on their GO MF annotations. See, for instance: https://www.ncbi.nlm.nih.gov/pubmed/26054620

Thats odd.  As a neuroscienctist I was totally unaware of this 'moonlighting'.  Are there gene ontologies that define protiens according to more fluid boundries?  Such as a hierarchical model?

Not sure what you mean. GO is a hierarchical DAG. Proteins with an annotation, also inherit all the parents of that annotation. In any case, my point was that it probably doesn't make sense to try or expect to be able to normalize across molecular functions. I wouldn't expect such classes to e very homogeneous in anything apart from function.

Forgive my ignorance, I'm hoping to enter a bioinformatics program.  My background is neuroscience I don't have a strong maths background.


text:
How can I find all the columns available in UniProt.ws package for R?
I am trying to find the subcellular localization of my 10,000 proteins using UniProt.ws package for R. However, I am unable to find all the columns available for query.
I used another package named Bioconductor which had easy way to lookup all the columns available using COLUMNS(org.Hs.eg.db)
I tried reading the documentation and couldn't find anything which can output subcellular localization. I am new to this field, any help is appreciated.

From the UniProt.ws manual:

columns shows which kinds of data can be returned for the
  UniProt.ws object.

So, you will need to go over the output of columns(your_uniprotws_object) if it includes the data type you are looking for.
And here is a related question.

Thank you. It worked for me now with this code:
up <- UniProt.ws(taxId=9606)
keys <- c("4957","81551")
columns <- c("PDB","GENES","SUBCELLULAR-LOCATIONS")
select(up, keys, columns, kt)

Do you have a small sample of the proteins you are looking for and what you have tried so far in R using Uniprot.ws ? (including output of it), it could be helpful to identify your problem and how to solve it.
Personally, I would rather consider Bioconductor as a platform for providing R packages dedicated to bioinformatics than a package to analyze data (but it's my opinion, maybe I'm wrong).

I have 10000 proteins which I am looking for sub-cellular localization. I tried to find the columns available in UniProt.ws using COLUMNS(UniProt.ws) but it didn't work. Finally, I looked into Bioconductor as you said and found the answer. https://support.bioconductor.org/p/78380/ . res <- select(up, "P35579", "SUBCELLULAR-LOCATIONS", "UNIPROTKB")


text:
Is there a Python package to convert InChi to molecular structures?
I am looking for a python package that can convert InChi keys (e.g.SGNXVBOIDPPRJJ-UHFFFAOYSA-N) to molecular structure or SMILES strings. Can I do that with RDKit?
The key comes from: https://www.metabolomicsworkbench.org/databases/refmet/refmet_latest.txt

Yes, RDKit can be used, however, if you installed it with conda it will not work out of the box for inChi key fetching. You can either spend some time installing the missing bit or use something else. It is not like Mol2 support in RDKit, which shouldn't be used, it is just not worth the bother.
I would suggest PubChem's RESTFul API.
Say
 import requests
 ethanol = 'LFQSCWFLJHTTHZ-UHFFFAOYSA-N'
 r = requests.get(f'https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/inchikey/{ethanol}/property/CanonicalSMILES/JSON').json()
 smiles = r['PropertyTable']['Properties'][0]['CanonicalSMILES']

Hey, thanks for the answer. Though I get {'Fault': {'Code': 'PUGREST.BadRequest',
'Message': 'Missing structure string',
'Details': ["An argument for 'inchi' must be provided in the URL path or POST data"]}} when I run your code.

Sorry for that. I fixed the typo by using your code. But I am glad it works now!

Slightly modified Matteos answer:
import requests
def get_smiles_from_inchikey(inchikey):
    r = requests.get(f'https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/inchikey/{inchikey}/property/CanonicalSMILES/JSON').json()
    return r['PropertyTable']['Properties'][0]['CanonicalSMILES']

inchikey = 'SGNXVBOIDPPRJJ-UHFFFAOYSA-N'
get_smiles_from_inchikey(inchikey)


text:
How I can change the name of multiple files at once in R or terminal?
I have 200 .vcf files in a folder with long names like
LP6005409-DNA_E03_vs_LP6005408-DNA_E03.snp.pass.vcf

How can I change the name of each file to, for example A1.vcf ....A200.vcf?
Renaming these files manually takes ages, help please.
> file.rename(list.files(pattern = "vcf"), paste0("A",1:length(list.files(pattern = ".vcf")),".vcf"))
Error in file.rename(list.files(pattern = "vcf"), paste0("A", 1:length(list.files(pattern = ".vcf")),  :
  'from' and 'to' are of different lengths

[fi1d18@cyan01 new]$ for i in *.vcf; do echo mv $i A${k}.vcf ; k=`echo "1+${k}" | bc`; done
mv LP2000325-DNA_A01_vs_LP2000317-DNA_A01.snp.pass.vcf A.vcf
(standard_in) 2: syntax error
mv LP6008031-DNA_A04_vs_LP6008032-DNA_D03.snp.pass.vcf A.vcf
(standard_in) 2: syntax error
[fi1d18@cyan01 new]$ for i in *.vcf; do mv $i A${k}.vcf ; k=`echo "1+${k}" | bc`; done
(standard_in) 2: syntax error
mv: try to overwrite `A.vcf', overriding mode 0500 (r-x------)?

just to complete your question, you can do it in R in one line.
After setting the directory as your working directory, just type the following:
file.rename(list.files(pattern = "vcf"), paste0("A",1:length(list.files(pattern = ".vcf")),".vcf"))

Maybe you can try file.rename(list.files(pattern = ".vcf"), paste0("A",1:length(list.files(pattern = ".vcf")),".vcf")). Also, if there is still an error, can you post the output of list.files(pattern = ".vcf") just to see the structure of your files in folder (maybe one file has a weird name not detected by the pattern).

Use perl-rename. This is usually called rename on Debian-based systems like Ubuntu or Mint, and perl-rename or prename on others. Assuming you have it as rename, you can simply do:
rename -n 's/.*\.vcf/"A" . ++$c . ".vcf"/e' *snp.pass.vcf

The -n causes rename to only print what it would do, without doing anything. So once you make sure it does what you want, run the command again without the -n.

On a Mac, you may need to install it first. Use brew install rename or use zmv in zsh (not bash) shell: https://superuser.com/q/152627/

@TimurShtatland note that there are many renames in the Linux world, just to keep things interesting ;)

If your filenames do not contain blank spaces, you can do it with a for-loop in bash:
Make a test run with echo so that the command is only printed but not yet executed:
for i in *.vcf; do echo mv $i A${k}.vcf ; let k++ ; done

And if the result looks good remove the echo statement to execute the mv commands:
for i in *.vcf; do mv $i A${k}.vcf ; let k++ ; done

Sorry I tried your code and I got error mentioned in the main post

@Angel note that there is no reason to use bc even if you do choose to do it this way. ((k++)) or k=$((k+1)) or let k++ all work to increment the value of a shell variable: k=0; while [[ $k -lt 10 ]]; do echo $k; ((k++)); done. Alternatively: for((k=0; k<=10;k++)); do echo $k; done.

Thanks for the comment!

EDIT:
(Based on the comments by @terdon, with minor changes)
Try this Perl one-liner. It is recommended for general use, even if your file names contain "unexpected" characters such as newlines.
perl -e 'for ( 0..$#ARGV ) { rename $ARGV[$_] => "A@{[ $_ + 1 ]}.vcf"; } ' *pass.vcf

Here, -e tells perl to look for code on the command line instead of a file with the script. @ARGV is the array of the command line arguments, here, files with old names.
for ( 0..$#ARGV ) { ... } loop goes through the indexes of elements of this array, assigning each implicitly to $_. The first argument to rename is the old file name: $ARGV[$_]. Instead of the comma separator of the method arguments I use => synonym for clarity in this method call.
The new file name is the second argument. @{[ ... ]} inside the double quotes is the Perl idiom to evaluate the expression inside, which is $_ + 1, to make the file numbering start with 1. If you don't care about this, it is easier to make it start with 0: just use "A$_.vcf", or, better, "A${_}.vcf". The last form needs to be used for cases like this: "A${_}pass.vcf" to prevent Perl from treating $_pass as a variable, where $_ is intended.
DEPRECATED:
Use the one-liner below when you are sure that the old file names do not have "unexpected" characters such as newlines. It is not for general use, see: Why you shouldn't parse the output of ls.
It will work with the example shown by the OP. It also works with blanks and tabs in the file names, but it breaks when file names have newlines.
The (deprecated) Perl one-liner:
ls *.pass.vcf | perl -lne 'rename $_, "A$..vcf";'

Here, ls | lists all files with old names, one per line and passes the stdout into stdin of the next command, which is the Perl one line script. -l removes the new lines from the input. -n reads the input one line at a time. These are all commonly used command line Perl options.
In the one liner script, $_ is the line read from stdin, and $. is the line number.

That is a slick perl one liner


text:
Should I use log2-CPM values (voom-limma) as input for my model?
We have created a model to integrate several OMICs data, but we realized that the maximum TPM values of RNA-Seq data were so big that had unexpected effects on our results. We hypothesized that this might be due to the negative binomial distribution inherent to this OMIC's layer.
To tackle this circumstance, we proposed to introduce as input the values obtained after using voom() (paper) in our count data. The reason behind that decision was that voom would transform the input data to be more accurately introduced to a linear model. There are, though, several issues/arguments against doing so:

The results obtained through voom() contain negative values, which do not make sense in our model. I tried to correct that by adding to all the results from a specific sample the absolute minimum value of that sample, but I am not confident whether that is a good way forward.
Original 0 values (non-expressed) turn into different negative values across samples (-6 in Sample 1, -4 in Sample 2), creating a difference that did not exist to begin with.
Some sources do not recommend to log-transform count data [1], including data from negative binomial distributions [2].

I am not really versed into distribution transformations, so any comment would be appreciated.

What sort of integration are you after? 2. Have you tried just using normalized counts?

The integration is supposed to give a better picture of transcriptomics as a whole, which would finally lead to specify the final destiny of a transcript based on its post-transcriptional regulation (whether it will translate or not). Till now we tried using TPM values to account for the different sequencing depths, as we have not used it yet to compare across conditions, but its behaviour across time. (Sorry, I erroneously flagged your comment while trying to upvote it :P)


text:
How to loop multiple function in shell script?
I need to extract sequences one after another consecutively from a large fasta files (multiple fasta files) and each extracted files to be saved in new fasta file (I mean the first sequence extracted from multiple fasta file to be saved as gene1 and second sequence extracted from multiple fasta file to be saved in another file and so on). For that, I tried a script like this,
for file in *.fasta do cat "$file" \
  | awk '/^>/ {printf("%s%s\t",(N>0?"\n":""),$0);N++;next;} \
          {printf("%s",$0);} \
     END {printf("\n");}' \
  | awk '{if ( NR==1 ) print $0;}' \
  | awk '{print $1"\n"$2}' | uniq > "new_${file}"

But, its not serving my purpose. Firstly, I could not loop my sequence extraction procedure from multiple fasta file and save as gene1 and gene2 and so on.
For example,

I have fasta files as follows,
org1.fasta
>locus1
ATGCGTAGAG
>LG12
TATGCATGAT
>bpl
TAGATGGAT

org2.fasta
>ysr
ATGTAGCGA
>inc
TGATCTATG
>pym
ATGCATAGA

org3.fasta
>siv
TAGTAGTAT
>tdp
ATGAGTCGA
>tem
ATGCTHATG

I need to extract first positioned sequences from all the org1.fasta org2.fasta and org3.fasta files and save in a new file gene1.fasta and second positioned sequences as gene2.fastaand third positioned sequences to be saved as gene3.fasta.
Expected output,

gene1.fasta
>locus1
ATGCGTAGAG
>ysr
ATGTAGCGA
>siv
TAGTAGTAT

gene2.fasta
>LG12
TATGCATGAT
>inc
TGATCTATG
>tdp
ATGAGTCGA

gene3.fasta
>bpl
TAGATGGAT
>pym
ATGCATAGA
>tem
ATGCTHATG

Please help me to do the same.
Thanks in advance.

You don't need a loop. You can do the whole thing with a simple awk one-liner:
awk -vRS='>' 'FNR>1{ printf ">%s",$0 > "gene"FNR-1".fasta"}' org*fasta

I ran this on your files and got:
$ for file in gene*; do echo "=== File: $file ==="; cat $file;  done
=== File: gene1.fasta ===
>locus1
ATGCGTAGAG
>ysr
ATGTAGCGA
>siv
TAGTAGTAT
=== File: gene2.fasta ===
>LG12
TATGCATGAT
>inc
TGATCTATG
>tdp
ATGAGTCGA
=== File: gene3.fasta ===
>bpl
TAGATGGAT
>pym
ATGCATAGA
>tem
ATGCTHATG

Explanation

-vRS='>' : set the record separator to >. This will make awk consider every fasta sequence as a single "line".
FNR>1{foo} : FNR is the line number of the current file. Because > is the first character of the entry, the first "line" will always be empty, so it needs t obe skipped.
printf ">%s",$0 > "gene"FNR-1".fasta" : this does the work. It prints a > followed by the current "line" ($0, this is the entire sequence + the header) and redirects it into a file called gene, then the current line number -1 (-1 is needed because of the issue mentioned above) and a .fasta extension.

=====
Note that org*fasta will be expanded with the default collation order. This means that org10.fasta will be processed before org1.fasta:
$ printf "org1.fasta\norg10.fasta\n" | sort
org10.fasta
org1.fasta

To avoid this and have the files in numerical order, you can use brace expansion instead of the *. So, if you know your files are all orgN.fasta where N can take values from 1 to 20, do this:
awk -vRS='>' 'FNR>1{ printf ">%s",$0 > "gene"FNR-1".fasta"}' org{1..20}.fasta

@K.Dineshkumar you're welcome! If one of the answers here solved your issue, please take a moment to accept it by clicking on the checkmark on the left. That will mark the question as answered and is the way that thanks are conveyed on the Stack Exchange sites.

Thank you so much @terdon. However, I am facing a problem, this script perfectly fine for up to 9 number of fasta files. If it exceed 9, it takes 10th fasta file and printing first instead of 1st fasta file sequence. Same way, it prints 20th fasta file second instead of 2nd fasta file. I do not know how to solve this.

@K.Dineshkumar ah yes, that is because of the order of collation. When you do echo org*fasta, it will sort alphabetically, so org10.fasta comes before org1.fasta (try with printf "org1.fasta\norg10.fasta\n" | sort). Try the updated answer.

Many thanks @terdon for your detailed explanation. With your explanation, I learned how to play around the script, in order to serve our purpose.

The title of your question indicates a shell solution, however you have added R and Python as tags. Here is an R solution making use of Biostrings pacakge (using two input files for simplicity):
library(Biostrings)

# create Biostring objects from each multifasta file
org1 <- readDNAStringSet("org1.fasta")
org2 <- readDNAStringSet("org2.fasta")

# loop over the Biostring objects
for(i in 1:length(org1)){

  output_file <- paste0("gene_", i, ".fasta")

  # create a Biostring object for each index
  genes_at_index_i <- c(org1[i], org2[i])
  assign(gene_name, genes_at_index_i)

  # save Biostring object as a multifasta file
  writeXStringSet(genes_at_index_i, output_file)
}

Output, gene 1:
>locus1
ATGCGTAGAG
>ysr
ATGTAGCGA

Output, gene 2:
>LG12
TATGCATGAT
>inc
TGATCTATG

Warning: The code above won't work if you don't have the same number of fasta entries in each multifasta file.

Thank you @haci for your help. However, I could not use this approach for 100 fasta files or above. Creating variable for each fasta file is tedious job.

You can do so in a loop, i.e. creating a list of Biostring objects with for or lapply and then looping over that list.

Thank you @haci, I would try to do the same.


text:
WGCNA module preservation error code with plotting summary statistics
I am trying module preservation with WGCNA using my own data. I am plotting the summary statistics including Zsummary, medianRank, observed IGP and IGP permutation p-value for two comparison data sets.
When plotting the scatterplots for IGP & p-value vs. module size I get following error code. I have re-run code and done practice data with no problems but unsure how to correct for this error? Does anyone have any experience with this and know how to resolve? Thank you kindly
 plotData[[2]][, 3] = cr[[ref]][[test]]$Acutal.IGP
  plotData[[2]],[, 4] = -log(cr[[ref]][[test]]$p.value + 1e -4)

errror in plotData[[2]][, 3] = cr[[ref]][[test]]$Acutal.IGP
number of items to replace is not a multiple of replacement length

Hi do you have the tutorial for plotting this?

My first guess is that there is a typo in your code, and Acutal.IGP should be Actual.IGP. I checked my tutorials and the word Acutal is in none of them, so I would assume you mistyped the variable name.

Oh my gosh! I am so sorry, I have re-read this so many times! and re-typed. I apologise I will get another pair of eyes to check as well as me before posting again. Thank you so much!


text:
AutoDock 4 -soft docking
How could I turn on/off the smooth potential in AutoDock4? I need to do two docking. The first one is soft docking with smooth potential and the second is re-docking without the smooth potential.

The smooth parameter is an optional keyword and command that changes the smoothing potential. If not provided, the default is 0.5 angstroms by default. Try passing 0 with this keyword in the AutoGrid Parameter File to turn off smoothing potential.
See the AutoDock 4.2.5 documentation for more details.


text:
Peak-calling using homer
I have a total of 78 ATAC seq samples from which Im trying to do a peak call .I tried this batchParallel.pl findPeaks peaks function but i couldn't find the output what it makes.
So I tried using loop but again the peak file is not written due to incorrect argument usage I guess.
These are my sample directories where i have all the tags
SRR2920564  SRR2920565

I used homer to make tag directories and now i using it for peakcalling.
here is my command.
for d in */ ; do findPeaks $d -style factor -o "${d}.txt" ; done

What i need is all the peakfile should be written as the directory name.
Such as if directory is "SRR2920564" then my peak file should be written in the same directory as SRR2920564.txt
I do run the loop it runs successfully but no output written. May be a mundane mistake am not able to figure it out.
Any suggestion or help would be highly appreciated .

Your output file is a bit weird.
findPeaks directoryname/ -style factor -o directoryname/.txt
You would be attempting to create a file called ".txt" inside each directory? (Also do you need those quotes).
Did you try running this with "-o auto" and checking in each directory if the corresponding files (i.e peaks.txt, regions.txt, transcripts.txt, etc..) were generated?

I tried running auto it gives an output peak.txt but when it would be easy if i can name them as the directory name for as unique name for each file

You just need to strip the / from $d, for which there are a number of options:
for d in */ ; do
    dname=`basename $d`
    findPeaks $d -style factor -o ${dname}.txt
done

or
for d in */ ; do
    dname=${d%%/}
    findPeaks $d -style factor -o ${dname}.txt
done

There are likely other ways one could go about this. Having said that, since you're not running these in parallel you could simply move the default output file (peak.txt) to a file name of your choice (using one of the strategies above).


text:
CRISPR/Cas9 screen analysis with Mageck: paired-end sequencing
I want to analyse data from a CRISPR/Cas9 screen (control vs. treatment) and I'm using Mageck (https://sourceforge.net/projects/mageck/). The problem is that I'm working with paired-end sequencing data, do I have to analyse R1 and R2 fastq files separately or is possible to do it at same time with Mageck?

The other answer is incorrect. The original poster is asking about paired-end sequencing data, not paired treatment/control experimental design.
Here is how you analyze paired-end sequencing data (this was an answer I gave on Biostars last year):
https://www.biostars.org/p/303656/


text:
How I move these files based on their names in R or terminal

EDITED

I have list of my samples in first column and corresponding file name for my .vcf files in second column of a file name clin_name.txt like below. For example name of one vcf is LP6008031-DNA_A04_vs_LP6008032-DNA_D03.snp.pass.vcf
OC/AH/120/E/T3  LP6005409-DNA_F01
OC/AH/126/E/T1  LP2000325-DNA_A01
OC/AH/136/E/T3  LP6005409-DNA_E02

This is the link of my file https://www.dropbox.com/s/mvzf47feiqe7x69/clin_name.txt.txt?dl=0
In the second file name pre_post.list I have information of pre or post treatment for each sample like this
pre-treated post-treated
OC/AH/120   OC/AH/065
OC/AH/126   OC/AH/069
OC/AH/136   OC/AH/074
OC/AH/160   OC/AH/084

This is link of this https://www.dropbox.com/s/ikrl1h8yoc0sngx/pre_post.list.xlsx?dl=0
As you suggested I typed
cat pre_post.list | awk '{print $2}' | sed 's@.*/@@'  > post_list.txt
cat pre_post.list | awk '{print $1}' | sed 's@.*/@@'  > pre_list.txt

Then I have
pre-treated
120
126
136
160
182

post-treated
065
069
074
084

in I typed ./test.sh in which I wrote
for item in  $(cat  pre_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /temp/hgig/fi1d18/TRG45/snp/
    else echo post_list $item;  fi;  done;

That returns
cp: cannot stat `*OC/AH/069/S/T3*': No such file or directory
[fi1d18@cyan01 snp]$

Or I run
for item in  $(cat  post_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /temp/hgig/fi1d18/TRG45/snp/pre/
    else cp *$out* /temp/hgig/fi1d18/TRG45/snp/post/;  fi;  done;

That returns
cp: cannot stat `OC/ST/026/S/T1*': No such file or directory
[fi1d18@cyan01 snp]$

I really don't know how to solve this :(

This is more of a job for snakemake or similar. You'd end up following the same logic in R if you really prefer that. Show us what you've tried so far, since there's no reason someone should write a custom program for you.

Assuming the names of all the files from one group are in file named list.
cat pre_post.list | awk '{print $2}' | sed 's@.*/@@'  > post_list.txt
cat pre_post.list | awk '{print $1}' | sed 's@.*/@@'  > pre_list.txt

Now you have two number list
    for item in  $(cat  pre_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /dir/for_pre_files/
    else echo post_list $item;  fi;  done

for post
    for item in  $(cat  post_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /dir/for_post_file/
    else echo Prelist $item;  fi;  done

or you could do this if one that isn't in post belongs in pre
    for item in  $(cat  post_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /dir/for_post_file/
    else cp *$out* /dir/for_pre_file/;  fi;  done

Sorry I tried your kindly code but did not work and I got this error in very  early stage  [fi1d18@cyan01 ~]$ cat pre_post.list | awk '{print $2}' | sed 's@.*/@@'  > post_list.txt
cat: pre_post.list: No such file or directory

You have to change the names in my bash code to reflect the names of the actual files on my computer

Thank you I am getting this the for loop  grep: clin_name.txt: No such file or directory
post_list 069 this being repeated for all of samples

You need to replace all the file names in the script with your file names.

Sorry what is "clin_name.txt" here? I mean where clin_name.txt comes from?  If I know what clin_name.txt is, I might be able to adopt that to my case Thank you

file that contains all OC/AH/120/E/T3  LP6005409-DNA_F01
OC/AH/126/E/T1  LP2000325-DNA_A01
OC/AH/136/E/T3  LP6005409-DNA_E02

Hello @Michael Hearn I tried what you kindly wrote for me edited in my main post but I am getting errors like cp: missing destination file operand after `/temp/hgig/fi1d18/TRG45/snp/pre/'


text:
In silico PCR: any models that can quantitatively predict the initial content from the end amplified result?
Let's say I've run the PCR and I now have a certain amount of the amplicon product. Is there any model that infers how much of it was there before the PCR started? I suspect it should be a probabilistic model (a distribution) parameters (shape) of which will strongly depend on the details of the PCR protocol and the characteristics of the amplicon...
If not, then may be there are forward models that predict the amount of the amplified result from the initial content?

You need to know the primer efficiency for this to work, do you know that?

@DevonRyan Assume, yes. You can assume I know everything. At this point I want to simply familiarize myself with any such model if there are any.

Have you considered just using a qPCR calculator such as this http://nebiocalculator.neb.com./#!/qPCRGen ? The model of qPCR is pretty simple really, so there are a few of these premade.

This is qPCR based on CT number, it is common place to identify patient viraemia for example

The basic equation for PCR kinetics (Equation 1) states that the amount of amplicon after c cycles (Nc) is the starting concentration of the amplicon (N0) times the amplification efficiency (E) to the power c. The PCR efficiency in this equation is a number between 1 and 2 (2 indicates 100% efficiency).  NCBI
In this  article  on research gate they said  amplifx  was the software you were looking for to, calculate "(i.e. expected efficiency in PCR reaction)". The amplifx even mentions "Testing primers in in silico PCR".