text:
missed cleavages combinations generation
I am looking for an R based solution to the following problem:
Trypsin does not cleave perfectly in proteomics experiments. So starting from a protein sequence FASTA I want to generate all possible peptides which COULD be generated if trypsin was allowed n number of "missed cleavages"
In informatics terms I would describe this as skipping a split during a strsplit() function
Example to clarify:
We start out with the protein sequence: GRGKA which is usually split by strsplit() (digested by trypsin) after every K and R
"GRGKA" %>% strsplit(., "(?<=K|R)", perl = TRUE)

[[1]]
[1] "GR" "GK" "A"

now imagine that trypsin did not cleave perfectly, for example missing always one K or R
"GRGK" "GKA"

would be generated in this case. Or if trypsin digested this protein really badly and would end up missing two cleavages we would end up with the sequence:
"GRGKA"

I would love to have a generalizable solution that allows to generate all combinations of sequences that arise with n number of missed cleavages.
If you need any more information please feel free to ask, sorry if I missed something!
jay.sf's requested example:
"FLGERTCTRNERQP"

perfect digest:
"FLGER" "TCTR" "NER" "QP"

one missed cleavage all combinations:
"FLGERTCTR" "TCTRNER" "NERQP"

two missed cleavages all combinations:
"FLGERTCTRNER" "QP" "FLGER" "TCTRNERQP"

three missed cleavages (there are a total of three K and Rs so trypsin would miss all)
"FLGERTCTRNERQP"

For longer sequences I think this becomes a frameshift like problem: If we count the K and Rs in the sequence:
FLGERTCTRNERQPFLGERTCTRNERQP
FLGE1TCT2NE3QPFLGE4TCT5NE6QP
FLGER_TCTR_NER_QPFLGER_TCTR_NER_QP (perfect digest)
FLGERTCTR_NERQPFLGER_TCTRNER_QP (one missed cleavage and always the uneven R and Ks)
FLGER_TCTRNER_QPFLGERTCTR_NERQP  (one missed cleavage and always the even R and Ks)
So I think for 1 missed cleavage we have two frames over the whole sequence and for two missed cleavages we would have three frames and so on....

Say "FLGERTCTRNERQP" would be a candidate, it has three possible splits. Can you show all possible results for 1, 2, and 3 skipped splits?

FLGERTCTRNERQP
perfect digest:
FLGER TCTR NER QP
one missed cleavage all combinations:
FLGERTCTR TCTRNER NERQP
two missed cleavages all combinations:
FLGERTCTRNER QP FLGER TCTRNERQP
three missed cleavages (there are a total of three K and Rs so trypsin would miss all)
FLGERTCTRNERQP

@KlemensFröhlich the answer below looks cool, you've enough rep to vote - definitely recommend it.

Essentially we want to combine 1, 2, ..., n - 1 splits, and if the variance of split lengths is greater than one we have a second solution. Here is a function that,  after the strsplit, first creates indices cb to later paste the respective splits together. I'm not quite sure what the maximum number of splits is, but this might be a start.
mccg <- function(x) {
  ## missed cleavages combinations generation ++++++++++++++++++++++++++++++++++
  gs <- \(n, k, length=n*k, right=FALSE) {
    ## helper, similar to `gl()` but w/o factorization
    if (right) {
      cbr <- gs(n, k, n, right=FALSE)
      rev(max(cbr) - cbr + 1L)
    } else {
      rep_len(rep.int(seq_len(n), rep.int(k, n)), length)
    }
  }
  ## split string after [K, R]
  spl_x <- el(strsplit(x, '(?<=K|R)', perl=TRUE))
  len <- length(spl_x)
  res <- if (len <= 1L) {
    ## if [K, R] at the end or missing
    x
  } else {
    ## nested list of combinations
    cb <- lapply(1:(len - 1L) |> {\(.) setNames(., . - 1L)}(), \(i, .len=len) {
      cb <- list(gs(.len, i, .len))
      if (i != 1L && var(tabulate(el(cb)))) {
        ## second solution, starting from right
        c(cb, list(gs(.len, i, .len, right=TRUE)))
      } else {
        cb
      }
    }) |> c(setNames(list(rep_len(1L, len)), len - 1L))
    ## paste according to sequences
    rapply(cb, \(x) vapply(unname(split(spl_x, x)), paste, collapse='',
                           character(1)),
           how='list')
  }
  return(res)
}

Usage
Two splits:
mccg(x='GRGKA')
# $`0`
# [1] "GR" "GK" "A"
#
# $`1`
# $`1`[[1]]
# [1] "GRGK" "A"
#
# $`1`[[2]]
# [1] "GR"  "GKA"
#
#
# $`2`
# [1] "GRGKA"

Three splits:
mccg(x='FLGERTCTRNERQP')
# $`0`
# [1] "FLGER" "TCTR"  "NER"   "QP"
#
# $`1`
# [1] "FLGERTCTR" "NERQP"
#
# $`2`
# $`2`[[1]]
# [1] "FLGERTCTRNER" "QP"
#
# $`2`[[2]]
# [1] "FLGER"     "TCTRNERQP"
#
#
# $`3`
# [1] "FLGERTCTRNERQP"

More splits:
mccg(x='FLGERTCTRNERQPFLGERTCTRNERQP')
# $`0`
# [1] "FLGER"   "TCTR"    "NER"     "QPFLGER" "TCTR"    "NER"     "QP"
#
# $`1`
# $`1`[[1]]
# [1] "FLGERTCTR"  "NERQPFLGER" "TCTRNER"    "QP"
#
# $`1`[[2]]
# [1] "FLGER"       "TCTRNER"     "QPFLGERTCTR" "NERQP"
#
#
# $`2`
# $`2`[[1]]
# [1] "FLGERTCTRNER"   "QPFLGERTCTRNER" "QP"
#
# $`2`[[2]]
# [1] "FLGER"          "TCTRNERQPFLGER" "TCTRNERQP"
#
#
# $`3`
# $`3`[[1]]
# [1] "FLGERTCTRNERQPFLGER" "TCTRNERQP"
#
# $`3`[[2]]
# [1] "FLGERTCTRNER"     "QPFLGERTCTRNERQP"
#
#
# $`4`
# $`4`[[1]]
# [1] "FLGERTCTRNERQPFLGERTCTR" "NERQP"
#
# $`4`[[2]]
# [1] "FLGERTCTR"           "NERQPFLGERTCTRNERQP"
#
#
# $`5`
# $`5`[[1]]
# [1] "FLGERTCTRNERQPFLGERTCTRNER" "QP"
#
# $`5`[[2]]
# [1] "FLGER"                   "TCTRNERQPFLGERTCTRNERQP"
#
#
# $`6`
# [1] "FLGERTCTRNERQPFLGERTCTRNERQP"

Simulation
Here are some toy data,
set.seed(42)
n <- 1e3
min_chr <- 5; max_chr <- 50
strings <- replicate(n, {
  repeat {
    smp <- sample(setdiff(LETTERS, c("B", "Z", "J", "O", "U", "X")),
                  sample(min_chr:max_chr, 1L), replace=TRUE)
    if (any(grepl('K|R', smp))) break
  }
  paste(smp, collapse='')
}) |> setNames(nm=_)

over which the execution of the function can be tested in a loop.
r <- lapply(strings, \(x) tryCatch(mccg(x), error=\(e) NA))

hi jay
Thank you very much for looking into this! This already looks promising, but the main goal would not be to get a function that can go as deep as possible but rather that can be limited at any level that the user wants to choose.
So basically it would be ideal to have a function that allows to set the missed cleavages to 0, 1, 2, etc and always receive all the combinations up to this point. Do you think your function could be modified to control the depth of allowed missed cleavages? Best Klemens

1) gsubfn is like gsub except the replacement can be a proto object containing a fun proto method which is a function whose first argument receives the proto object and which inherits a count (and other) built-ins.  The match is input to the function as the second argument and is replaced by the output of the function.  Here we replace each R or K with itself and if it is more than the k-th occurrence then with a semicolon as well.  Then split by semicolon.
library(gsubfn)

# inputs
x <- c("GRGKA", "GRGKAK")  # vector of inputs
k <- 1 # no of splits to skip

p <- proto(fun = function(self, x) paste0(x, if (count > k) ";"))
strsplit(gsubfn("[RK]", p, x), ";")

giving
[[1]]
[1] "GRGK" "A"

[[2]]
[1] "GRGK" "AK"

2) Insert a semicolon after each R or K and then remove the first k.  strsplit what is left.  This uses base R only.
xx <- gsub("([RK])", "\\1;", x)
for(i in seq_len(k)) xx <- sub(";", "", xx)
strsplit(xx, ";")
# same output as above


text:
Using PANDA to infer gene regulatory networks without ppi
I'm using this tutorial (https://github.com/netZoo/netZooM/blob/master/tutorials/panda/panda.pdf) to use PANDA to infer gene regulatory networks, but without ppi values (which I've been told can be done easily using MATLAB). I replicated the codes given in the tutorial but I omitted all the portions with ppi. After doing so, I'm getting a lot of warnings and errors that makes me think I've made syntax errors in the process of taking out the ppi portions. Here's what I have:
addpath(genpath('path/to/netZooM'))

load("myMat.mat")
myTable = struct2table(myTable);
addpath(genpath('path/to/netZooM'));
fileparts(which('netZooM'));

exp_file = 'tests/test_data/expression.txt';
motif_file = 'tests/test_data/motifTest.txt';
panda_out = '';
lib_path = 'C:/Users/ual-laptop/Documents/bioinformatics.m/netZooM';
save_temp = '';

alpha = 0.1;
save_pairs = 0;

expressionData = readTable('expression.txt');
size(expressionData)

[TF, gene, weight] = textread(motif_file, '%s%s%f');
TFNames = unique(TF);

disp('Reading in expression data!');
fid = fopen(exp file, 'r');
headings = fget1(fid)
n = length(regexp(headings, '\t'));
frewind(fid);
Exp = textscan(fid, ['%s', repmat('%f', 1, n)], 'delimiter', '\t');
fclose(fid);
GeneNames = Exp{1};
Exp = cat(2, Exp{2:end});
[NumGenes, NumConditions] = size(Exp);
fprintf('%d genes and %d conditions!\n', NumGenes, NumConditions);
Exp = Exp';

disp('Computing coexpression network:');
GeneCoReg = Coexpression(Exp);

disp('Reading in motif data!');
[TF, gene, weight] = textread(motif_file, '%s%s%f');
TFNames = unique(TF);
NumTFs = length(TFNames);
[~,i] = ismember(TF, TFNames);
[~,j] = ismember(gene, GeneNames);
RegNet = zeros(NumTFs, NumGenes);
RegNet(sub2ind([NumTFs, NumGenes], i, j));
fprintf('%d TFs and %d edges!\n', NumTFs, length(weight));

disp('Normalizing Networks:');
RegNet = NormalizeNetwork(RegNet);
GeneCoReg = NormalizeNetwork(GeneCoReg);
TFCoop = NormalizeNetwork(TFCoop);

disp('Running PANDA algorithm:');
AgNet = PANDA(RegNet, GeneCoReg, TFCoop, alpha);

AgNet = panda_run(lib_path, exp_file, motif_file, panda_out, save_temp, alpha);

sum(AgNet, 2)
sum(Agnet, 1)

Are there any glaring syntax errors in my code? I'm thinking my errors could also stem from maybe cloning the github repository into the wrong place (which is a separate topic) but it's definitely possible that there are syntax errors since this is my first time doing this kind of stuff.

source: Stack Exchange

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

sort: votes


text:
How to do average macrophages signature boxplots, with Bulk RNA data from 2 control and 2 treatment samples and LM22 signature matrix from cibersort?
I have bulk data from 4 samples of triple negative breast cancer, which 2 of them are control and the other 2 are treated. The goal here is to estimate the average macrophages signatures, for M0,M1 and M2 in each of the 4 samples.
My bulk data:
Gene    MCtrl1_TPM  MCtrl2_TPM  MTreat1_TPM MTreat2_TPM
KHDC1L  0.0 75.15   0.0 0.0
CD24    1.53    41.56   0.11    0.06
MAMLD1  5.2 11.94   0.0 0.0
SMIM31  4.1 10.24   0.0 0.06
C16orf95    7.9 13.24   0.0 0.08
OCLN    2.69    7.06    0.06    0.03
RBM14-RBM4  16.97   1.95    0.0 0.0
RPS10-NUDT3 0.0 0.0 31.91   12.06
NDUFC2-KCTD14   15.06   49.98   0.41    0.0
ATP9A   2.7 10.05   0.01    0.0
UST 0.38    10.63   0.0 0.01
LOC100996709    12.25   30.58   0.0 0.0
ADAM8   10.13   2.14    0.03    0.0
ADGRG2  0.05    21.49   0.0 0.0
MBNL2   4.61    5.04    0.0 0.0
PIGBOS1 56.73   51.33   0.6 0.09

To do that, I thought was to use the CIBERSORTx software more specifically, the option - Impute Gene Expression, High-Resolution Mode (tutorial 5) with my bulk data and signature matrix file, the LM22 signature matrix file of 22 immune cells provided by CIBERSORTx, but with the selected columns correspondent to the 3 types of macrophages.
Signature matrix file:
Gene    M0  M1  M2
ABCB4   7.951360487 27.55771004 121.4322774
ABCB9   56.25188749 29.38087813 40.39026017
ACAP1   48.97221651 30.97425375 23.05515036
ACHE    13.6929629  1373.06962  48.58576071
ACP5    13644.1723  1062.475909 7633.960182
ADAM28  54.96648703 274.639997  360.6820442
ADAMDEC1    4464.262133 8524.022801 5004.969419
ADAMTS3 45.56049098 32.34016681 33.35142773
ADRB2   143.0178779 787.6239741 363.1850058
AIF1    1841.269508 1611.228571 11610.29035
AIM2    181.9571    2965.888162 148.2845223
ALOX15  6.621226619 5.412216318 430.3353048
ALOX5   788.5397751 83.30089783 248.3724724
AMPD1   76.54456664 62.39764428 59.95269372
ANGPT4  39.33257003 21.218127   11.8097012
ANKRD55 123.600429  96.20392388 48.55417359
APOBEC3A    85.2400852  4501.140479 148.7993403
APOBEC3G    278.7681406 1036.993881 258.3004281
APOL3   238.5980956 11009.29392 437.5272083
APOL6   13.10943681 967.6451456 7.784612131
AQP9    4583.500675 4309.341662 399.4746192
ARHGAP22    175.4574354 351.5922058 272.9890572
ARRB1   333.3962724 946.0138184 395.5562302
ASGR1   108.7218916 25.23626149 155.4913385
ASGR2   36.11978652 26.3738691  760.7257707
ATHL1   52.25820981 167.5032595 132.631945

But I encountered this error, trying to do that:
## Formatted for easier reading
Error in runCIBERSORTxGEP(weights, mixture, classes, QN, maxsamples, label, :
Abort! The number of mixture samples is not greater than the number of cell types.
Calls: CIBERSORTxHiRes -> CIBERSORTxGEP -> runCIBERSORTxGEP Execution halted

How this error is possible, when I have 3 types of macrophages and 4 samples? There's another way to estimate the average of the 3 types of macrophages signatures? Specifically in a graphic like this?If yes, how?

Welcome to the site. The raw data has 4 columns (2 controls) the processed data has 3 columns. Why is that? Specifically what do M1, M2 and M3 mean? Is this macrophage1, marcophage2, ... 3 . How does MCtrl1_TPM  MCtrl2_TPM  MTreat1_TPM MTreat2_TPM relate to M1, M2 and M3?

Thank you! This image presented is an example from another article that i read and wanted to do that type of graphic to my analysis. Yes M0,M1 and M2 are macrophage0, macrophage1, macrophage2. Given that i dont have single cell data of my data to exactly estimate average of each macrophage type in each sample, o do that i used a signature matrix from cibersortx software, LM22 signature matrix file of 22 immune cells, and from that matrix i selected only the columns correpondent to the macrophages types, that i'm interest to work, and i then tried the Impute Gene Expression, Group Mode.

Okay I see. This is bioinformatics but if you don't get a response, you could try Biology SE because it's formal package that for example an immunologist would know how to use. Just mentioned you've tried here first.

I have a question: CIBERSORT and CIBERSORTx give me same results? thnk

You can use IOBR R package for this purpose. IOBR is an R package to perform comprehensive analysis of tumor microenvironment and signatures for immuno-oncology.
if (!requireNamespace("remotes", quietly = TRUE))
install.packages("remotes")

if (!requireNamespace("IOBR", quietly = TRUE))
remotes::install_github("IOBR/IOBR")

library(IOBR)
library(tidyverse)

file = read.delim(file.choose(),row.names=1)        #bulk data from 4 samples

# estimating the average macrophages signatures using cibersort method
cibersort<-deconvo_tme(eset = file, method = "cibersort", arrays = TRUE, perm = 100 )

head(cibersort)

#barplot
res<-cell_bar_plot(input = cibersort[,c(1,15:17)], title = "CIBERSORT Cell Fraction")

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


text:
Is there a data analysis software (free) or code that I can use to view normalized CDR3 amino acid length distributions?
Is there a data analysis software (free) or code that I can use to view normalized CDR3 amino acid length distributions?
Regarding code, preferably using Python or R.
EDIT: Added a picture of the non-normalized CDR3 amino acid length distribution. The issue is that ES0017_metastasis has more samples than all other individuals which is why it has the highest peak at every amino acid length from 10-20. Was wondering if normalization was possible? Plot generated via Immunarch package in R.
EDIT2: Now able to access OriginPro and Graphpad Prism via university.
EDIT3:

EDIT4:

Here are the results, what do you think...
The first chart is the log transformation (one extra as well, don't know how to delete). The bottom chart is the min-max transformation.
The x-axis refers to amino acid length. Each discrete block is an amino acid length, e.g. 12 amino acids long, etc.

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Welcome to the site. It is unlikely free software is available, although PAML did capture this distribution, I think for HIV (thats about trees). The question is not clear as to what the goal of normalisation is. A histogram would be useful. Normalisation is usually requires trial and error based on the histogram. It is likely to be over dispersed but understanding the distribution is a good start.

Apologies for the lack of clarity, edited the quesiton.

Final update
Okay the answer below is based on the second graph from the top. Looking at the update the correct answer is closer to the first answer presented. EDIT 3 is normalised data. It appears that the y-axis is probably natural log transformed (loge) or log2 and this transformation has normalised it.
Thus there is not need to normalise the data, the program has already done it in EDIT3 and in the very top graph now presented.

This looks impressive, thanks. I get it. There are two way to do this. If you are normalising the frequency (y-axis) is what you are wanting to normalise then either (in this instance its the wrong answer [see below the ----]):

a log10 transformation (even Excel will log this data)
a max-min transformation (thats in Python, likely in R - I dunno a general package maybe SPSS??, dunno)

You then test via the normal distribution using with Shapiro–Wilk test or Kolmogorov–Smirnov test
I suspect the max-min transformation will work. There are other ways in but it will involve discarding data that is >3SD from the mean. Try the above and see how it works.
Note, a specialist statistician will request you perform a QQ plot, because the above tests of normalisation are not perfect.
This is not standardisation this is conformation to the normal distribution.

This could be about bin size if you have a discrete range for gene length and doesn't need to be very big - normally you need to try a few values, in this case its just 1 amino acid bin ... This can be used as a sliding window. It will pull the data into a normal distribution. That may not be very well explained, but I'm pretty sure that will work.
1 amino acid sounds weird but part the failure to see a normal distribution is a result of your calculation: you are splitting an amino acid into fractions and thats whats causing the irregularities in the data. Thus between values 10 and 20 on the x-axis (presumably size range difference) you should only have 20 bars - but you've got loads > 100 which isn't possible.
Do that try SW and KS tests and you're sorted.


text:
Constructing a pangenome graph using multiple samples
I am working with yeast strains (specifically S. cerevisiae) that have been cultured and evolved over many generations in the lab. I am using the reference S288C strain when I do certain analyses but I want to be able to construct my own reference using the "in-house" strains. Ultimately I want to be able to call variants when mapping strains to this "in-house" reference to make it easier to visualize/understand the difference in strains.
I am trying to use VG-tools to create a pangenome graph using 10 different strains all sharing the same reference and I'm honestly a little confused on how to go about doing it. I was able to make a "vg" file for each sample individually but I know downstream I want to be able to map other strains to this pangenome as a reference and call variants to it rather than the standard reference sequence.
I'm mainly looking to see if anyone has built a "pangenome graph" using multiple samples and how they went about doing it. If anyone could give me some tips/tricks/advice for this kind of task that would be great!
I have a phylogenetic tree of all the strains but how would I use that to see the differences in variation?

Hi @rimo what do you mean by a "graph". In graph theory a "graph" is a network. If its a classical bivariate plot, what is the y-axis? Finally, what about simply doing a tree?

Could you post a tree - remove the taxa labels to ensure confidentiality - and I will tell you how to read it. The answer is in the tree. Where the SNPs are located is a different analysis BTW.

@M__ rimo is talking about a genome graph, like the example shown in the vg repository

@gringer I did raise this point originally. This is graph theory, but graph theory is based on trees (it provides the direction).


text:
Can the PGS Calculator give the contribution of each variant to the overall PGS?
I am using the PGS Calculator (https://github.com/PGScatalog/pgsc_calc/) to compute the polygenic scores of individuals based on a scoring file.
PGS Calculator works wonderfully and outputs the aggregated score for each individual.
However, I am interested in the contribution of each variant to the score of each individual.
At the moment, I have:

IID
PGS

33432984732
10

32432809751
20

I would like the following table:

IID
Variant 1
Variant 2
Variant 3

33432984732
3
4
3

32432809751
6
8
6

in which I can see the contribution of each variant to the overall score.
Is this possible with PGS Calculator? If not, is there a tool other than PGS Calculator that I can use? Or should I use my own scripts?

This would be a hard feature to support natively, because some PRSs have 10 million+ SNPs.
There is a tool that does this, but I don't want to recommend it because it will require debugging their source code upon running, and it is less accurate than what you are using here.
I suggest doing --scorefile "my_custom_scores/*.txt" if you are interested in only a few variants. Add only one score to this file. This may cause some error especially if some of your samples do not have all positions defined.
If you have more variants, there is a way to combine multiple scores in a single scorefile instead of a directory.


text:
`bedtools`: Intersecting/map into only one interval when there are multiple matches?
I have the following list of intervals

Chr
start
end

chr1
0
5

chr1
5
10

chr1
10
15

chr1
15
20

chr1
20
25

chr1
25
30

and the following CpGs I'd like to intersect with the intervals and aggregate the last two columns using sum.

chr
start
end
end -1
coverage

chr1
1
3
2
1

chr1
4
6
5
1

chr1
7
9
8
2

chr1
10
12
11
2

chr1
13
15
14
1

chr1
16
18
17
2

chr1
19
21
20
2

chr1
22
24
23
2

chr1
25
27
26
1

I ran bedtools map -a cgis.bed -b s.bed -c 4 -o sum, but when a CpG falls
across two regions, bedtools counts it in both. Here the 5 - 10 intervals should have only the second row of the CpG table: a sum of 2. However, since the second row spans 5, the 4 - 6 interval falls into both the 0 - 5 and the 5 - 10.

chr
start
end
bedtools sum coverage
correct sum coverage

chr1
0
5
2
2

chr1
5
10
3
2

chr1
10
15
3
3

chr1
15
20
4
4

chr1
20
25
4
2

chr1
25
30
1
1

Can I use bedtools to place it in the lower interval only and not both? Intersects should be start inclusive, but end exclusive. This does work as I'd like if I reduce every end column element in the second CpG table by one.

What question are you trying to answer exactly? The second row of your CpG table really does overlap with two rows of your first table.

@bricoletc I want to pick the first of the overlapping regions, not both. If there are two the second one should not be counted. I think its a question of end-inclusive vs not. If every end in the CpGs table is reduced by one, I get the correct answer. I was wondering if its possible to get this correct answer without editing the input CpG file.


text:
What is the required input peak file format for IDR (Irreproducible Discovery Rate)
I have managed to install IDR (https://github.com/nboley/idr) correctly (idr --version returns 2.0.3), but I can't get it to run on my peak files (*.bed) without crashing.
My peak files are generated using HOMER (findPeaks) and they are converted to proper *.bed format via HOMER/pos2bed.pl.  But when I run idr on my *.bed peak files as described in the IDR guide page like below:
idr --samples peak1.bed peak2.bed
it always returns an error inside the function bed_load() :
signal = float(data[signal_index])
IndexError: list index out of range
I suspect it's related to peak files (likely columns) not being in the format the program idr.py is expecting, and I couldn't find much about the required "input" file format in the idr help pages.
Can someone provide examples of an idr.py friendly peak file, maybe just a couple of rows with the proper/expected column types and numbers? Thank you in advance.
I am including the guesstimate guidance from ChatGPT, which is not necessarily correct, but maybe a pointer in the right direction - it's suggesting the IDR format has to be consistent with narrowPeak formats having 10 columns, and the actual column that IDR uses in it's analysis is column 7 - which it calls it the "signal value". But what exactly is this measure (the peak strength, the peak height etc.)? Any insight would be appreciated? Also, on where this type of information may have been discussed, seeing how critical it is to running the tool successfully.

The IDR readme on usage describes very clearly what the OUTPUT format is, but is not clear what the INPUT format is (column descriptions) - which is what one needs to successfully form the INPUT and run IDR. My typical HOMER/pos2bed.pl converted row is less than 8 columns which at this point, I am certain is the cause of IDR failures. I am attempting to figure out which column IDR actually uses, and which are not essential for it's considerations (they just get duplicated in the output).
Here are the error messages I get:
https://github.com/nboley/idr/issues/41

Please show the first few lines of your input files, and (if possible) the full command you've run, and the entirety of the error message.

As mentioned in the IDR readme on usage, there are examples in the IDR github repository here. Here's one of them:
chr16   31085416    31086009    .   0   .   362.546677850757    -1  4.12368834166759    250
chr12   109490129   109490582   .   0   .   360.329496957805    -1  4.12368834166759    215
chr13   41345146    41345561    .   0   .   358.624215585062    -1  4.12368834166759    207
chr19   50083722    50084001    .   0   .   266.603941224671    -1  4.5777559727315 186

From the looks of it, it seems that the output is a line-filtered version of the input format (i.e. output and input formats are identical).
I have no idea how IDR works (never used it before), but the error you've added to your post suggests there's a range mismatch. My guess based on that error message is that the two input files don't have overlapping specified regions.


text:
Errors in AutoDock
I am using Autodock and I'm getting an error in the step where we get the glg file (i.e. Run--> Autogrid) I am not sure what should I do...Please help me if someone knows how to tackle this. Plus yep there's no space in the folder file name and I have all exe and dat files in the same one folder.
                         $Revision: 1.109 $

Compilation parameters:  NUM_RECEPTOR_TYPES=20 NEINT=2048

MAX_MAPS=16 NDIEL=16384 MAX_ATOM_TYPES=14

energy_lookup table has   655360 entries of size 8

e_vdW_Hb table has   401408 entries of size 8

Maximum number of maps that can be computed = 16 (defined by MAXMAPS in "autocomm.h").

GPF> npts 18 12 24                        # num.grid points in xyz

Number of grid points in x-direction:    19

Number of grid points in y-direction:    13

Number of grid points in z-direction:    25

GPF> gridfld 98 Alanine_over_Phenylalanine.maps.fld # grid_data_file

C:/Bioinformatics/Docking/autogrid4.exe: ERROR:  Grid data file needs the extension ".fld" for AVS input

C:/Bioinformatics/Docking/autogrid4.exe: Unsuccessful Completion.

Welcome to the site @Khanishka. One thing thats not clear is whether the installation is successful. Have you successfully run test data?

@M__  Yes installation is successful. As far as I know yes I have


text:
Adding BAF and Log2 ratio to VCF file
I have access to pre-processed VCF files from a genetic dataset, and am trying to use the MoChA tool to calculate mosaic chromosomal alterations in this dataset. However, MoChA requires BAF and log2ratio values for each of the variants, which are missing from the VCF files given, and the original genome array CEL files are unavailable so I cannot redo the variant calling and include these statistics
Instead, they have provided us with nexus format CNV files (from the Axiom CNV summary tools) for each of the samples with a host of log2 ratio and BAF files across the different Axiom probes.
I am just a bit confused about how we can use these files to import the BAF and log2ratio values into the VCF for MoChA. Any suggestions?
A couple sample rows of the CNV file are included here (I've removed the sample identifier and replaced with ---).
Sample  ProbeSetName    Chromosome      Position        Log2Ratio       BAF
---     AFFX-KIT-000001 9       101258881       -0.440535940556599      0.482727669963674
---     AFFX-KIT-000002 4       164934874       -0.106227164241716      0.057035393165642

Welcome to the site. The manual from GitHub https://github.com/freeseek/mocha appears to describe this in detail both for BAF and log2 (which the manual calls LRR for Log R Ratio). I would recommend reading the manual because it is extensive and looks very clear.


text:
stats.sh error in BBmap package
I am trying to calculate the N50 value from the assembled FASTA file. I used stats.sh from the BBmap package. I executed the following command
./stats.sh  in=/Home/Downloads/BBMap_39.01/bbmap/resources/contig.fasta

I got the following error.
Exception in thread "main" java.lang.RuntimeException: Input file does not appear to be valid: /Home/Downloads/BBMap_39.01/bbmap/resources/contig.fasta
    at jgi.AssemblyStats2.process(AssemblyStats2.java:248)
    at jgi.AssemblyStats2.main(AssemblyStats2.java:40)

How to resolve this error?
The sequences are shown below. My file also contains some lengthy sequences.
>NODE_1_length_3013_cov_4.540906
TGGGGGATCATGTAACTCGCCTTGATCGTTGGGAACCGGAGCTGAATGAAGCCATACCAA
ACGACGAGCGTGACACCACGATGCCTGTAGCAATGGCAACAACGTTGCGCAAACTATTAA
CTGGCGAACTACTTACTCTAGCTTCCCGGCAACAATTAATAGACTGGATGGAGGCGGATA
AAGTTGCAGGACCACTTCTGCGCTCGGCCCTTCCGGCTGGCTGGTTTATTGCTGATAAAT
CTGGAGCCGGTGAGCGTGGGTCTCGCGGTATCATTGCAGCACTGGGGCCAGATGGTAAGC
CCTCCCGTATCGTAGTTATCTACACGACGGGGAGTCAGGCAACTATGGATGAACGAAATA
GACAGATCGCTGAGATAGGTGCCTCACTGATTAAGCATTGGTAACTGTCAGACCAAGTTT
>NODE_2_length_2292_cov_16.637908
CAAAGGAATAATAATTGACTGATGATTGGTTCATTAATTACTTCATTAGTTAATAGATTG
CAGTCTTTATAAATAAGTTTATAGTTGTATTCCACAGAATTTCTGATCCTGTCTCTGCCC
TACCCAATGTTTTATCATTAAACTTCATTTACATAAGATAGTATGTTCAACAAATATAAA
GATTTGCTGGAATGGAAGGCACAGTAAGATATTTCACAGGATAAAATTAAAGTAGATGTC
AACAGAGTAAATCAAAGTGAAATTTAATAGAGTTAAACAGTAAGTGTTTTAAATACCAGG
ATTAAATTCCATGAAGCTGCAGCTTTTAAAAAAAGGCACAAAGTGGAATAAGAAAAAATG
ACTTTCAGAAAGTTAATGTAAAGAAAGCTCACGGTTTTCATTGGCTGCCTGCTCAGTTTG

Also, why you run this on a file provided by bbmap rather than your own?

@ATpoint Thank you for your comment. contig.fasta is my file. I just copied to the resources directory.

Hmmm. Is this a permissions error? Can you type ls -l /Home/Downloads/BBMap_39.01/bbmap/resources/contig.fasta and report the output?

@M_. Thanks for your comment. ls -l shows the list of files. The first two lines are shown. -rwxrwxrwx  1 user user  2674 Jun 12  2020 addadapters.sh
-rwxrwxrwx  1 user user  2674 Jun 12  2020 addadapters.sh

I agree with @atpoint as best guess. ls /home/Downloads/BBMap_39.01/bbmap/resources/contig.fasta , i.e. lower case h.  It could be other things. Can you try their test data?

./stats.sh in=notExistingFile.fa
Exception in thread "main" java.lang.RuntimeException: Input file does not appear to be valid: notExistingFile.fa
        at jgi.AssemblyStats2.process(AssemblyStats2.java:248)
        at jgi.AssemblyStats2.main(AssemblyStats2.java:40)

The error comes when the file does not exist. My best guess is that it is a typo on your side. Is the path really /Home with capital "H"? Most OS like Linux and macOS use lowercase. Anyway, make sure the file exists.


text:
obtaining a fragment data file bed from bam
Im trying to follow this pipeline https://github.com/epifluidlab/cragr on cfDNA WGS. I have bam files from paired end reads aligned with bwa mem.
The workflow requires bgzipped and indexed fragment data files as input.

I have been searching how can I obtain BED files from fragments coordinates instead of reads coordinates, and came up with no answer.
Does someone know how can I obtain such bed file?

Hi, what do you mean? :) I don't know how to do it. I'm familiar with bamtobed from bedtools but that gives me a bed file based on reads. I need a bed file based on fragments instead :)

I see so you are looking whether there is a bedtools option for fragment output. Got it, thanks. Would be good to modify the question with the info.

There's an example in the README (https://github.com/epifluidlab/cragr#example-datasets):

For you information, below are steps from the BAM file to the
fragment file:

samtools sort -n -o qsorted.bam -@ 4 sample.bam
samtools view -h -f 3 -F 3852 -G 48 --incl-flags 48 qsorted.bam |\
  bamToBed -bedpe -mate1 -i stdin |\
  awk -F'\t' -v OFS="\t" '{if ($1!=$4) next; if ($9=="+") {s=$2;e=$6} else {s=$5;e=$3} if (e>s) print $1,s,e,$7,$8,$9}' |\
  sort -k1,1V -k2,2n |\
  bgzip > sample.frag.bed.gz
tabix -p bed sample.frag.bed.gz


text:
running fragpipe
I am getting the following error when I run FragPipe on my Ubuntu machine 22.04 how do I fix this?

Thanks for your question. As a general point it's always good to write out the error in asci rather than have an image. It makes it difficult to read and remember at least on my browser. The context of the question is processing mass-spec and I don't know if anyone here knows (maybe someone does). Is there a file format issue per chance?

the file format that I use as input is control.raw and treated.raw for fragpipe and affter running it i got the following log which is below:

Checking database...
Failed in checking /home/ms/MS/Control_1.raw. Will ignore it.
Failed in checking /home/ms/MS/Treated_1.raw. Will ignore it.
Batmass-IO binaries for Thermo support and/or Thermo native libraries not found
Batmass-IO binaries for Thermo support and/or Thermo native libraries not found
There are corrupted files. Please remove those files and re-start the task.
Checking spectral files...
/home/ms/MS/Control_1.raw: Scans = 0
/home/ms/MS/Treated_1.raw: Scans = 0
Process 'MSFragger' finished, exit code: 1
Process returned non-zero exit code, stopping

Edit your post and add the plain text content there with proper formatting.

Thanks for your valuable time, well my error has been resolved through Fragpipe support Here I am adding the discussion link below:
https://github.com/Nesvilab/FragPipe/issues/1207


text:
Reference genome (bwa -mem)
I want to do mapping with bwa mem. I downloaded a reference genome from ncbi database which is (Danaus plexippus). I got one ncbi_database.zip file which contain 2 folders and 3 files of doc,assembly_report_json file. I read one file named read.md which came after downloading but its not helping. What should be done in this case?
I don't have an issue with using bwa-mem; the problem is with the reference genome from NCBI. I downloaded it from https://ncbi.nlm.nih.gov/datasets/genome/GCF_009731565.1 but after that I checked the downloaded file and I got 5 folders into it and I'm confused which one to use. 2 files have .fna extension and the rest are reports that are in the json file format.

Welcome to the site @moon. My role is to try to get the questions as clear as possible. Firstly it's worth noting bwa-mem2, but you want to try getting bwa mem working first. The read.md will just be the installation instructions. Are you looking for the manual? Is the following helpful, https://hcc.unl.edu/docs/applications/app_specific/bioinformatics_tools/alignment_tools/bwa/running_bwa_commands/ and follow the bwa mem instructions. Maybe try it and then post back on specifically what the problem is, e.g. the error message?

The answer by @atpoint is the best answer to the question as the question is written right now. My advice is to consider accepting this and further questions can then be asked. This Lepidoptera genome assembly and a reference genome is definitely needed.

Only accept answers if they solve the problem you're having. If the answers don't match your problem, feel free to clarify the question to explain further.

I downloaded genome from https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_009731565.1/ I got 5 files in it so if you are confused which one to use in that so GCF.fna having genomic data with fna extension is used as an reference genome if you are also new to this things this might be useful.
Thank you

Please accept the other answer as well, as it addressed your question and was given before you added your own answer and accepted just the one you added

What you need is the genome sequence in fasta format, typically suffixed .fa, for example from https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_009731565.1/
This file needs to be indexed via bwa index.


text:
what does 8 stands for in bwa mem command line
What does this command implies in bwa mem
bwa mem -t 8 index/referencegenome.fna read1.fastq.gz read2.fastq.gz > output.sam"""

What does -t 8 stands for here?
Thank you

-t is an important parameter here and is important to understand. bwa is a CPU calculation and the more CPUs the job (calculation) is spread across the faster the calculation will take to complete and the less chance of the calculation exceeding the RAM capacity of the machine. "Exceeding RAM capacity" (I don't know the technical term), is when the process "grinds to a halt" and never completes.
Basically what can happen is the indexing step instead of taking ~3 hours, could easily be still running after 24 hours if the machine has little RAM and -t isn't set sufficiently high.
One CPU is known as "one core".
For a regular desktop if you set -t X where X is number of CPUs - 1 in my opinion that is good practice. So if the machine has 8 CPUs then X = 7. If you are working on a cluster, e.g. GCP or AWS, then X is the maximum number of cores (CPUs) available.

thank you .you actually got my question this really did helped me rather then just knowing that it just stands for threads

If I'm using a regular desktop I should write it as Bwa mem -t 1 rather than 8 am i right?

pedantic note: the indexing in bwa is not multithreaded, only the alignment (bwa mem) is.
you should set it as high as your machine allows to finish in time. note that the higher "-t" goes, the more memory the process will consume. there is no fixed rule for this since in bwa the memory usage depends on the reads you're aligning, its complexity and length, but usually something in the range from 1-16 with available memory somewhere in the range of 8-32gb works fine. if memory is limited then keep "-t" down, something like "-t 4" will probably do.

thanks I'm gonna prefer 8 for now then

-t 8  is good. It is good to test empirically, so if it takes >3 hours indexing occurs try a different value. I do not recommend -t 1 at all. The point is that -t is a pivotal parameter for RAM usage. There is a difference of opinion on multithreading versus RAM. I agree that multithreading followed by single threading is common (I've written code that does that) and could be an issue (it's really complicated).

When you have questions on software then the first thing is to do your due diligence and read the manual:
https://bio-bwa.sourceforge.net/bwa.shtml
Then simply checking for the "-t" argument will tell you that it is the number of threads used for the alignment process.

sorry but I was not asking about what is t ...I was asking why to put 8 in that which is not given in the manual I wanted to know if I can change it or is there a specific rule about it and yes I'm definitely new to this bioinformatics stuff but I'm not new to science I'll learn things anyways thanks for replying

Then you should've asked "how to set -t" and which implications does choice of "-t" has.

my bad i was not aware of that after -t it's the thread number


text:
low mapping percentage in bwa mem
After aligning paired-end reads to a reference genome, I am getting low percentage:
34636166 + 0 in total (QC-passed reads + QC-failed reads)
34346076 + 0 primary
0 + 0 secondary
290090 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
5706057 + 0 mapped (16.47% : N/A)
5415967 + 0 primary mapped (15.77% : N/A)
34346076 + 0 paired in sequencing
17173038 + 0 read1
17173038 + 0 read2
3807784 + 0 properly paired (11.09% : N/A)
3996748 + 0 with itself and mate mapped
1419219 + 0 singletons (4.13% : N/A)
177284 + 0 with mate mapped to a different chr
65388 + 0 with mate mapped to a different chr (mapQ>=5)

What has to be done to increase the mapping percentage?
I'm attaching the fastqc report before trimming and after trimming too
first-  R1 without trimming
second- R2 without trimming
third- forward paired
fourth- reverse paried
rest are overrepresentations
After using trimmomatic the quality is getting more bad and even the mapping is of low percentage
Mash screen result in the last image

My initial guess would be that you have a low quality reference genome, or you mixed up the reference genome, or you have significant contamination of your library, or that your sequencing run is very low quality. I'd suggest running fastqc and posting the results here, posting some statistics about your genome, and generally explaining where your reads are coming from and what your reference is.

@MaximilianPress Thanks for your reply . I have the data of  Parantica aglea I done fastqc that report was showing probelm of gc content and base quality even after trimming the sequence reads with trimmomatic it got more worse

I am unclear if this solved your problem, possibly consider posting the fastqc outputs if you want more help. Seems like the reads might just be bad but hard to evaluate from what you say.

@MaximilianPress I edited the question added the fastqc report before and after trimming.

The reads honestly look fine to me from a sequence quality standpoint. FASTQC thinks that the basecalls are good based on the green light. Your first few nucleotides are a little weird compositionally but that's probably due to UMIs or barcodes in your library, without seeing more information. A read mapper would just soft-clip those. We are not seeing e.g. an excess of Gs that you would expect from an instrument failure. You are not showing the overrepresented sequences plot which is probably most likely to show any systemic issues, but none of the plots that you are posting look bad.

@MaximilianPress so what exactly to be done like BWA MEM is showing the low percentage how to improve it.

I'd suggest going through the hypotheses in my first comment. The reads look at least sort of ok (though you should check the overrepresented sequences tab of Fastqc). So you could look at 1) contamination, try using mash screen 2) wrong/bad reference genome (you will know better how to do this, maybe mash screen the reference genome against the screen database?) 3) sample swaps of your library. I don't really have more to say than that.

@MaximilianPress thank you I'll follow as I alreadly used different softwares for trimming and mapping but results are same no changes like fastp,bowtie2,bwa-mem2 so maybe mash screen will help to know it better I edited my question I added up overrepresentation screen of fastqc take a look on it .

@MaximilianPress I did mash screen against my sample and this was the last output. I added the snap in the question will you take a look at it

Just answering, following extensive back and forth in the comments.
I think that your reads are fine:

FASTQC doesn't see anything too weird. You see some noise at the first few basecalls but I am not too alarmed about that.
FASTQC finds overrepresented only a telomeric repeat and its reverse complement (TTAGG). That is something you'd expect to find in WGS reads.
Mash screen finds that the only match, with a small but non-trivial number of k-mer matches, is the monarch butterfly Danaus plexippus, e.g. another lepidoptera. So that makes sense.

Therefore, the issue is likely that your reference genome has issues or has been swapped out with something else.
Where did your reference assembly come from?

Did it get mixed up somehow?
What does it look like when you run QUAST?
What does it look like when you run mash dist of your assembly against the same mash database? Is it maybe a bunch of endosymbionts or other contaminants?
What happens if you map your reads against a high-quality lepidopteran assembly like the monarch butterfly or this?
What happens if you map your assembly against a high-quality lepidopteran assembly? Is it syntenic?

I got the reference genome from ncbi database. I tried out mapping with another reference genome of same subfamily which is Danaus chrysippus and the result are same and there are no reference genome available apart from this so I'm thinking might be there is some genetic variation between the parantica aglea and Danaus plexippus . So, de novo assembly will the solution for it ? or is there any another solutions which I can look for ?

@Moon It would be helpful if you could specify the reference assembly initially in the future. I still don't know which reference assembly you are working from, but I am guessing that it is not in the same genus as your reads. Yes, you can expect significant variation even between fairly closely related organisms- 10-20% is not crazy unless they are pretty close. Yes, I would do a de novo assembly. Even with just Illumina you are likely to do ok compared to these other references.

I have the sample sequencing reads of Parantica aglea and the reference genome I'm using is Danaus plexippus which comes from the same subfamily I was trying to do Reference-guided genome sequence assembly pipeline so I initially started with fastqc and then trimming and here I'm on bwa mem where there I'm getting low percentage even after changing the reference genome and software tools the story is still same as you said might be the issue with reference genome so de novo assembly can be done I thought. Might this will clear what I'm working on

@Moon I would suggest not doing a reference guided assembly. I would instead do a de novo assembly. I expect that the issue is that your genome and the Danaus genome are substantially different in structure and you will not be able to find anything meaningful by reference guided assembly, and may be actually misled. After you perform the de novo assembly you can still try to scaffold the resulting contigs using the Danaus reference using a tool like RagTag.


text:
Why doesn't bedtools getfasta pick up the first base?
I am using bedtools getfasta to convert genomic coordinates to sequences. I feed it a bedfile containing say chr1 s.start s.end. When I check the sequence output, it does not include the s.start base.
For example if the bedfile entry is chr1 23456 23678, it gets the sequence from 23457-23678. Is this a feature and I should use a tag to make it include the first base?
The command I am using is bedtools getfasta -fi genome.fa -fo gene.fa -bed input.bed
I tried this in bedtools versions 2.27.1 and 2.30.1

BED files are 0-based. I link this old but gold biostars thread for more background. Essentially, from normal counting with start 1 you subtract 1 from the start coordinate.
For example in the sequence ATCGGG the TCGGG substring would have "normal" (aka human-readable) coordinates 2-6, but in BED format that would be:
chrSomething    1    6

Here is the bedtools "proof":
$ cat foo.bed
chrSomething    1       6

$ cat foo.fa
>chrSomething
ATCGGG

$ bedtools getfasta -fi foo.fa -bed foo.bed
>chrSomething:1-6
TCGGG


text:
Renumbering atoms in PDB using OpenBabel
I'm trying to add polar hydrogens to protein (PDB ID 3l3l cleaned from water and ligands) and convert it to PDBQT using OpenBabel by:
obabel -ipdb target.pdb -opdbqt -O target.pdbqt -xr --partialcharge gasteiger --addpolarh
It's working well, but all hydrogens are appending to the end of file. I want to renumber atoms so that residues in PDB would go in increasing order and hydrogens were located "between" corresponding heavy atoms, like AutoDock Tools do this.
Is there a way to renumber them using OpenBabel? Or any other command-line utilities?

Try reopen in biopythown and save pdb  with preserve_atom_numbering = False

Please add a small bit of the pdb with non ordered H so that people can try to help you without having to spend too much time

See https://www.mail-archive.com/openbabel-discuss@lists.sourceforge.net/msg05790.html

I see the prob openbabel ( at least this one https://www.cheminfo.org/Chemistry/Cheminformatics/FormatConverter/index.html  OPEN BABEL 3.1.1) doesnt name the Hs correctly (their are all just H) so that biopython parser drops them as duplicates. Use pymol h_add to a pdb then save directly from pymol or run conversion in openbabel

mmmmh even my pymol assigns wrong names to H atoms compared to the ones in the https://www.ebi.ac.uk/pdbe-srv/pdbechem/ or http://ligand-expo.rcsb.org/ld-search.html

https://phenix-online.org/documentation/reference/ready_set.html ReadySet! is a program designed to prepare a PDB file for refinement as in ReadySet! Refine!!!. It will add hydrogens to the protein model using phenix.reduce and to the ligands using eLBOW. The appropriate restraints are also written to disk. Hydrogens can also be added to water molecules. Deuterium atoms can be added to facilitate dual xray-neutron refinement. Metal coordination files are also generated.

https://kpwulab.com/2021/02/19/pymolchimeraadd-hydrogens-to-your-structure/ chimera should work fine


text:
Discrepancy with featurecounts analysis using a forward stranded and reverse stranded protocol
My RNAseq analysis pipeline is as follows: fastqc (read quality is good, some overrepresentation of adaptor sequence) → trimmomatic (trimmed adaptor sequence, qc report after trimming suggests the overrepresented adaptors are gone) → HISAT2 (used a forward stranded alignment protocol) → featurecounts
My HISAT 2 command line:
hisat2 --rna-strandness F -x HISAT2/grch38/genome -U "$fastq_file" | samtools sort -o "HISAT2/${sample_name}_aligned.bam"

My featurecounts command line:
featureCounts -s 2 -a hg38/gencode.v44.chr_patch_hapl_scaff.annotation.gtf -o "quant/${sample_name}_featureCounts.txt" "$bam_file"

However, the strange thing is that if I use a forward stranded protocol in feature counts (with the forward stranded alignment I generated with HISAT2), the assignment rate is 5%, with the majority of unassigned reads being having no features. But if I use a reverse stranded protocol in feature counts, the assignment rate is 70%. Does this suggest my data is reverse stranded?
My library was prepared by the NEB directional RNA library prep kit that claims to produce forward stranded libraries, so I’m really confused by this result. How should I proceed with this data? Should I run an alignment with a reverse strand protocol and just proceed with a reverse strand assignment? Or should I proceed with the forward alignment and forward strand assignment data because that’s what the library is supposed to be?
Thank you so much for your help!

Not a direct answer but an alternative approach: The nfcore/rnaseq nextflow pipeline infers strandedness automatically and feeds in the "right" parameters to the tools used in the pipeline so you do not have to worry about this. You can even find the inferred strand info int he output directory of the nextflow run. And if you want to check your reads, you can give a try with Salmon or RSeqQC.

This kit uses the dUTP method to preserve strandedness. That (to my knowledge) always produces a reverse library. Not sure where you have your information of a forward library from, but we used this kit for many years and it is reverse.


text:
seeing full results on R
running this code on R, using the Biostrings package,
vmatchPattern("GAAGACCGAGGCCACGCGGAGTACGAACGAGGGTACAGTCTTC", mySeq)

in the console, it writes:
...
<22 more elements>

How can I see all the results? The print function returns the same result.

It returns an IRanges object. You can transform it to a data.frame using as.data.frame() and then either write as a file with something like write.table() to view outside of R or use View() in RStudio to explore it.

Thanks a lot :)

Accept the answer to provide closure to the port. @FarzadBeikpour


text:
How to create a phylogenetic tree from diverse mitochondrial genomes
I would like to create a phylogenetic tree for the most species in my dataset.
I'm starting with around 1200 species, but since it's not good practice to align short and long sequences I tried filtering only for the species with the [15k-18k bp] range in length, but this leaves me with around 350 species in total.
Those sequences vary a lot in size, here's an histogram of the number of species I have data for

As you can see I have a lot of sequences with length [12k-15k bp] and [15k-18k bp] (values on the x axis are scaled by a factor of 3000 and then rounded).
I would like maximise the species genetic information.
To do so I would like to extract the part of genome that is shared the most among all the species.
Is this achievable by using a regex?
I am not sure on the length and type of the query to be sure of not extracting non-homologous regions as it would create a big problem later on for the phylogenetic tree creation.
I can't make an alignment for all species because of the lack of computational power. That is why I'm trying to subset a piece of DNA that is the most shared (prior to the alignment, but I don't know if that is doable or logic)
My question is: how would you proceed to get the most phylogenetic data out of this dataset?
If you would use a regex, what type of query would you use?
More Info
Alignment
When I try to align 10 sequences drawn from the [15k-18k bp] range I get this result (the numbers on the right are the TaxID of the species).

We can't see much here so I used an online MSA viewer
[https://alignmentviewer.org/]

Data
My data is one Fasta file containing mitochondrial genomes for all the species.
To do this preliminary analysis I only used 10 randomly chosen species with a sequence of length in the range [15k-18k bp].

Have you considered using only genes for the alignment instead of the entire sequence?

Is there a way for you to find an existing vertebrate phylogeny (this seems to be your group of interest) based on a larger set of sequences and to prune it to the set you care about? e.g. idk here, see .newick and .tre dowloads. iTOL makes interactive tree exploration/pruning pretty easy.

Your overall approach is cool, given caveats (below) to adapt the method. I can fairly authoritatively say your specific methodological approach will not work regarding its phylogenetic strategy without changing it. This is due to a phenomenon where the number of mutations becomes so great it is impossible to recover the ancestral state, basically it randomises the phylogenetic signal, this especially occurs at the 3rd codon nucleotide positions, i.e. 1 in every 3 nucleotides, 33.3% of the entire data set. This violates a central tenant of phylogenetic theory, therefore I am afraid that you cannot build a tree using this approach, but again you can adapt it to work (below).
If you want only to use nucleotides then 16S is the only target available. However, 16S has been a good target for amphibians, therefore just aligning that gene and building a tree is cool. Its perfectly solid given prior work in vertebrates ...

Chan, K.O., Hertwig, S.T., Neokleous, D.N. et al. Widely used, short 16S rRNA mitochondrial gene fragments yield poor and erratic results in phylogenetic estimation and species delimitation of amphibians. BMC Ecol Evo 22, 37 (2022).

If you wish to use protein coding genes you need to switch amino acids and align those. In terms of conservation there alot of established tools rather than "reinvent the wheel". One program you could consider is AAcon, available here. This would give you can empirical basis for deducing which gene is

the most conserved
the most represented, i.e. "read-depth", or species representation.

... after its aligned, thereby yielding a genetic marker for phylogenetics. It uses a large number of metrics therefore most scenarios would be accommodated.
COX1 is good target and likely gene to emerge from AAcon analysis.
Building protein phylogenetic trees is a separate question.

Lacking computational power to align. That is solvable, even on small computers, and you can place that as a separate question stating your computer capacity (RAM and CPUs). It's not difficult to align 1000s even on a small machine, it's simply the strategy to do this. Algorithm performance and accuracy change rapidly.

It looks like these are random mitochondrial sequences pulled from NCBI, which is probably why the alignments aren't working well.
As a starting point, I recommend you make sure that all the mitochondrial DNA sequences are rotated and reverse-complemented as necessary to make sure that the same gene matches position 1 of the sequence. I prefer doing this for COX1, because it's a well-conserved gene. After doing that, it'll be much easier to look at the alignment of that gene alone, which will give you an idea of the complexity of extending the task to the full length of the mitochondrial genome.
I don't expect you'll find a long exact-match sequence that is consistently conserved across all species, even within COX1. It wouldn't surprise me if the longest conserved sequence were less than 16bp, with a few repeats of that sequence throughout the genome. If that is the case, you need to think about what to do with that information. Given the complexity of this intermediate problem you've asked about, it would be a good idea to think first about what you would do after you discovered there were no such sequence.


text:
how to install cibersort in R, to do deconvolution using my RNA-seq Data with public single cell data
I have RNA-seq Data from lung cancer immunotherapy patients but i don't have any single cell data from my sample, so i wanted to do a deconvolution with a public single cell data. I thought of doing the analysis with the cibersort website plataform but i wanted to obtain the results in R, to be easy to work with them.So i followed the steps in this article, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5895181/ :
library(e1071)
library(parallel)
library(preprocessCore)
source("CIBERSORT.R")
results <- CIBERSORT(‘sig_matrix_file.txt’, ‘mixture_file.txt’, perm=100, QN=TRUE)

But in when i did source("CIBERSORT.R"), i obtain this:
"Error in file(filename, "r", encoding = encoding) :
  cannot open the connection
In addition: Warning message:
In file(filename, "r", encoding = encoding) :
  cannot open file 'CIBERSORT.R': No such file or directory"

I don't understand i was supposed to have a file?

Is the R script available on the Cibersort website? I was looking for it some time ago but couldn't find. I don't have my credentials for the site on this computer so can't check now.


text:
error in downloading the annotation of sequences on NCBI
I have developed an R function to download the "organism" name of the numerous NCBI accession numbers which are already saved in an object.
The function:
get_organism_names <- function(accession_numbers) {
  organism_names <- sapply(accession_numbers, function(accession) {
    record <- entrez_summary(db = "nucleotide", id = accession)
    return(record$organism)
  }, USE.NAMES = FALSE)
  names(organism_names) <- accession_numbers
  return(organism_names)
}

Although it works for some accession numbers, for others it encounters an error stating for example:
Warning: Esummary includes error message: Invalid uid KU253616.1 at position= 0Error: No esummary records found in file

How can I either solve this error or ask R to skip the problematic accession and keep downloading the rest?

Have you tried running it in a loop with a single accession each time instead of giving it all of them?

Also, please give us all of the code. Presumably, you are loading some libraries, which ones? We need example code that we can use to reproduce the problem. See https://stackoverflow.com/help/minimal-reproducible-example. We need example input, and then code that we can run locally that reproduces the error.

You could use the tryCatch function to handle inputs that throw an error. A good example on how to write them can be found here. Though the message indicates it's a warning, so does the function stop running when you get this message?


text:
GO Enrichment Analysis
I want to annotate 25,000+ proteins with GO_ID, TERM, and CATEGORY prior to doing enrichment analysis. Can I extract this information by submitting my sequences? Meanwhile, my species is a fish, and it's a non-model organism. Since my species is a non-model species, I have annotated the proteins from the whole genome by myself using interproscan.  There I was able to find GO IDs but there was no term or category which is not helpful for Enrichment analysis.
Seq type:
>maker-ctg000210_np1212-snap-gene-104.0-mRNA-1 p/t AED:0.26 eAED:0.29 QI:38|0.66|0.5|1|0.33|0.25|4|0|342
MPCPSKRSRAAKKRKADLRSTEVPSKRSPVAQVDLVSCDEVGRPRKKAHTQSSDAIQSSDCQLQCQLQCQ
LQCQLQSSDAIQANSTNSQTNSTNSLNMCHDNNATDCCFLNESVDFPERSVLMGSFHQGDSRFGNVRNKQ
CGAISLTAVLKSKMKNVWSWETGDLDDVLIKGTVLYSTPRFEVTGVTAVTHATAQMTDSSSCATSSVNYT
VPGKDVKQNVCRPLYVHEVGVTLNRDCASKDENRSCECNVRKPETKQDKQSFNESTTDDVILVNVETHKQ
SLKRNLQATDEVIISDVRFHKQVFKPLTRHDQDALCTCLDLQNENVYLGLMATTDLDDMGCP

After annotation:
maker-ctg001000_np1212-snap-gene-144.133-mRNA-1       710fdf6795f6f6bfbb960a017c734103        1139    ProSiteProfiles PS50056 Tyrosine specific protein phosphatases domain profile.  354     408     12.21363        T       20-02-2023      IPR000387       Tyrosine-specific protein phosphatases domain   GO:0016311|GO:0016791

But I want to use my protein sequences and the following output I need:
Seq_id    GO_ID    GO_term    GO_categoty

I would appreciate any of your suggestions. Thanks.

Please [edit] your question and show us an example of your input data and the output you would want from those data. Are these known, annotated proteins? What IDs do you have? Are these just raw sequences with no annotation? Do you need to assign GO terms based on homology with other, known sequences?n

That looks like an unknown, novel protein. Does it actually have GO annotations? How are you annotating it?

It's unclear what you're asking. Can you please give some more information about what you will do with this output once you've found it.
The single annotation you have provided has GO IDs, and using AmiGO, I get the following results:

GO:0016311

Name: dephosphorylation
Ontology: biological_process

GO:0016791

Name: phosphatase activity
Ontology: molecular_function

Bear in mind that GO is a hierarchical annotation, and most GO enrichment analysis programs ignore this hierarchy when carrying out their analysis, which is even more important for a less studied non-model species, where annotation is unlikely to be complete at any level.
Using a standard enrichment program that ignores the hierarchy would require fully populating the GO tree for each protein. For example, using the inferred tree view from AmiGO for that one protein:

GO:0016311 -> GO:0008150|GO:0009987|GO:0008152|GO:0044237|GO:0006793|GO:0006796|GO:0016311
GO:0016791 ->  GO:0003674|GO:0003824|GO:0016787|GO:0016788|GO:0042578|GO:0016791

Consider this when thinking about the type of output you want. This will not be a single GO term/category per sequence; each sequence will have multiple GO IDs, multiple terms, and multiple categories, with many overlaps and lots of gaps.
Further information about exploring the GO database (both automated and manual) can be found here:
https://amigo.geneontology.org/amigo/software_list


text:
How do I select a subset of genes for functional enrichment(GO/KEGG) analysis from WGCNA results?
I am working on using WGCNA for a bulk RNA-seq experiment. I have three experimental conditions (control, treatment 1, and treatment 2) and have a total sample size of 18 (6 per group). Using the online WGCNA tutorials, I could analyze my datasets after variance stabilizing transformation (vst). I got the summary of network analysis results as a CSV with Gene significance (GS), p.GS, module membership(MM), and associated p-value (p.MM).
I want to extract a set of significant genes from the module with the highest positive correlation to my treatment conditions. As I am just getting acquainted with WGCNA interpretation, which of these parameters in the CSV should I use for filtering further (GS, p.GS, and module membership and associated p-value (p.MM))?
Also, how would I choose the limits of the criteria? Should it be using the scatter plot (MM vs GS)?
As I am pretty new to WGCNA, any help is appreciated. Thanks in advance!

Consider filtering on every parameter: Gene Significance (GS), p-value for Gene Significance (p.GS), Module Membership (MM), and p-value for Module Membership (p.MM).
This doesn't mean you need to necessarily filter on every parameter, but it depends on what exactly you want to do, and you might reasonably choose any/all.

For GS, you may want to ensure you're capturing genes strongly associated with your experimental conditions.

For p.GS, you may want to ensure statistical relevance, (so filter out genes with high p-values). Similar for p.MM.

For MM, you may want it to be high to make sure the genes you select are strongly co-expressed and are relevant members of the module.

Plot MM vs GS to see if you can see an obvious filter threshold.


text:
how to use BLAST+ in an optimal way
I have a 32-bp query that I want to perform BLAST+. my query represents different RNA viruses and I do expect that my hits get more than 500,000.
but in practice, when I perform it I get around only 700 hits. would you please kindly help me how to increase the efficacy of my BLAST+?
I leave my syntax here.
blastn -task megablast -db nr -outfmt 7 -remote -max_target_seqs 99999999 -perc_identity 80 -query D:/Query.fas  -out D:/result.txt

Can you explain why you expect many queries? More importantly, what is the biological question you are trying to answer here? Are you looking for distant homologs? Close homologs? Why did you choose to use megablast? That's an odd choice for such a short sequence.

Thanks for your comment.
I expect many hits because there are many hits for my query as I know my query is common among RNA viruses, it is such a conservative seq.
Yes, I want to extract all sequences in the NR database that are similar to my query.
I think megablst was a bad choice, because when I choose -task blastn my hits were significantly improved (+500,000), but still I encounter an error by blast+ stating "search cannot be formatted" resulting in finding just some of those hits.

Please [edit] your question and explain this. Megablast is only useful for very similar sequences. And long ones, usually. It doesn't really make sense to use it for a 32bp query, that will be smaller than the smallest thing it can find.  So please [edit] and explain what you are trying to do. You don't want similar sequences, presumably you are looking for homologs, right? I mean explain the biological question you are trying to answer, so we know how to help.

Megablast is designed to find very similar sequences. It is highly specific but has a very low sensitivity. This is because it has a very large word size. By default, I think it will use a word size of 28 nt. That's almost your entire sequence, so if your target sequence doesn't have 28 contiguous nucleotides that are identical to 28 contiguous nucleotides of your query, no hit will be found.
If you are looking for distant homologous sequences, you should use simple blastn which defaults to a word size of 11. That will still likely be too small, you might want to change it to something lower (your blast tool should have an option for that), but it will give you more results than megablast.
Finally, if these 32nt are translated into protein, you can consider using tblastn instead.


text:
Error in snakemake when running conda envs in container
I am facing an issue when running snakemake with a container and conda.
Let me explain better.
I am developing a workflow (locally) in snakemake, where some rules use a conda env (defined in /envs/env.yaml) and other rules use wrappers (so they also are based on conda envs).
An example here:
rule clean_spike:
    input:
        sample_ref="{}results/bam/{{id}}.tmp.bam".format(outdir),
        sample_spike="{}results/bam_spike/{{id}}_spike.tmp.bam".format(outdir),
        sample_ref_index="{}results/bam/{{id}}.tmp.bam.bai".format(outdir),
        sample_spike_index="{}results/bam_spike/{{id}}_spike.tmp.bam.bai".format(outdir),
    output:
        sample_ref="{}results/bam/{{id}}.clean.bam".format(outdir),
        sample_spike="{}results/bam_spike/{{id}}_spike.clean.bam".format(outdir),
    conda:
        "../envs/pysam.yaml"
    log:
        "{}results/logs/spike/{{id}}.removeSpikeDups".format(outdir),
    script:
        "../scripts/remove_spikeDups.py"

Locally, the pipeline runs smoothly.
Now, I wanted to test the pipeline on another machine (server from the institute), and so I run the command
snakemake --containerize > Dockerfile
I got a dockerfile and  I generated an image that was uploaded on dockerhub here.
At this point, I copied the workflow on the other machine, I added to the Snakefile the line:
containerized: "docker://davidebrex/chiprxsnakemake:latest"
and tried to run the pipeline with: snakemake -c 10 --use-conda --use-singularity.
Ideally, snakemake should run the conda envs within the docker container for each rule
The problem is that i get an error related to the fact that conda is missing..
I checked the docker image and I can run conda and see the installed envs.
Any suggestions about what is going on?
Thank you!
This is the error I get:
/usr/bin/bash: conda: command not found
Traceback (most recent call last):
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/__init__.py", line 792, in snakemake
    success = workflow.execute(
              ^^^^^^^^^^^^^^^^^
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/workflow.py", line 1246, in execute
    raise e
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/workflow.py", line 1242, in execute
    success = self.scheduler.schedule()
              ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/scheduler.py", line 636, in schedule
    self.run(runjobs)
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/scheduler.py", line 685, in run
    executor.run_jobs(
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/executors/__init__.py", line 158, in run_jobs
    self.run(
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/executors/__init__.py", line 582, in run
    future = self.run_single_job(job)
             ^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/executors/__init__.py", line 635, in run_single_job
    self.cached_or_run, job, run_wrapper, *self.job_args_and_prepare(job)
                                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/executors/__init__.py", line 623, in job_args_and_prepare
    self.workflow.conda_base_path,
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/workflow.py", line 456, in conda_base_path
    return Conda().prefix_path
           ^^^^^^^
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/deployment/conda.py", line 653, in __init__
    shell.check_output(self._get_cmd("conda info --json"), text=True)
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/shell.py", line 61, in check_output
    return sp.check_output(cmd, shell=True, executable=executable, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/subprocess.py", line 466, in check_output
    return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/subprocess.py", line 571, in run
    raise CalledProcessError(retcode, process.args,
subprocess.CalledProcessError: Command 'conda info --json' returned non-zero exit status 127.

I tried to change the snakemake command and remove --use-conda but then I get an error related to the fact that the tools are missing.
I checked out other questions on stackoverflow but they were not helpful in my case (like this one: here).

I have zero experience with --containerize, but "I checked the docker image and I can run conda" versus "conda: command not found" when running the workflow really smells like a conda shell setup issue to to me, like people run into when trying to call conda in a script versus interactively.  (No idea where to point you from there, though!)

Thank you. At the end I changed snakemake version and it worked..quite weird

Actually now i still get the same error, I also opened an issue on snakemake github, I am clueless


text:
Calculating p-value and adjusted p-value for pre-normalized microarray data with fold change precalculated
I have a microarray dataset with two mutants dataset that has already been normalised, and the fold change values for each gene in each mutant versus the wild type have been calculated. I'm interested in determining the statistical significance of these fold changes.
What I've tried: I initially considered using limma, but I wasn't sure if it was appropriate for pre-calculated fold change values.
What I was expecting: I was hoping to obtain a list of genes with their corresponding p-values and adjusted p-values, indicating their significance. Ideally, this would help me identify genes that are significantly upregulated or downregulated in the mutant compared to the wild type.
Given that I only have the fold change values for two biological replicates (and no raw intensity or count data), what would be the best approach to achieve this? Any guidance, recommended tools, or R packages would be greatly appreciated!
Thank you in advance!

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.


text:
bcbio run not running samples in parallel
Sorry for the cross post but I'm not getting a response on bcbio github:
I'm trying to run bcbio on 20 samples and it's only running one at a time as it steps through the yaml seemingly. I'm using slurm. I was wondering if someone can look at this bcbio snippet and if they see anything wrong?
Slurm script
#!/bin/bash

#SBATCH --job-name=bcbiopipeline_GPC3_PDx
#SBATCH -t 1-00:00
#SBATCH -c 1
#SBATCH -p core -n 32
#SBATCH --mem-per-cpu=40G
#SBATCH --error=job.%J.err
#SBATCH --output=job.%J.out
#SBATCH --export=ALL

echo "Load bcbio..."
module use /projects/ngs_prodifx/local/software/modules/
module load bcbio-nextgen/1.2.9

bcbio_nextgen.py seqc.yaml -n 72 -t ipython -s slurm -q core -r t=0-72:00 -r conmem=20 --timeout 3000 --tag "rna"

yaml snippet:
resources:
  default:
    memory: 30G
    cores: 72
details:
- algorithm:
    aligner: star
    disambiguate: mm10
    expression_caller:
    - salmon
    - kallisto
    quality_format: standard
    strandedness: auto
    transcriptome_fasta: /swiftcache/ngs/oncology/analysis/todd_creasy/working/reference/GRCh38.primary_assembly.genome_GPC3CART.fa
    transcriptome_gtf: /swiftcache/ngs/oncology/analysis/todd_creasy/working/reference/gencode.v29.annotation_GPC3CART.gtf
  analysis: RNA-seq
  description: 10_LI6612_Baseline_8743
  files:
  - /swiftcache/ngs/oncology/analysis/todd_creasy/working/fastq/SM-745/SM-745-10_R1_001.fastq.gz
  - /swiftcache/ngs/oncology/analysis/todd_creasy/working/fastq/SM-745/SM-745-10_R2_001.fastq.gz
  genome_build: hg38
  metadata:
    panel: Baseline
- algorithm:
    aligner: star
    disambiguate: mm10
    expression_caller:
    - salmon
    - kallisto
    quality_format: standard
    strandedness: auto
    transcriptome_fasta: /swiftcache/ngs/oncology/analysis/todd_creasy/working/reference/GRCh38.primary_assembly.genome_GPC3CART.fa
    transcriptome_gtf: /swiftcache/ngs/oncology/analysis/todd_creasy/working/reference/gencode.v29.annotation_GPC3CART.gtf
  analysis: RNA-seq
  description: 11_LI6612_Baseline_8765
  files:
  - /swiftcache/ngs/oncology/analysis/todd_creasy/working/fastq/SM-745/SM-745-11_R1_001.fastq.gz
  - /swiftcache/ngs/oncology/analysis/todd_creasy/working/fastq/SM-745/SM-745-11_R2_001.fastq.gz
  genome_build: hg38
  metadata:
    panel: Baseline
- algorithm:
    aligner: star
    disambiguate: mm10
    expression_caller:
    - salmon
    - kallisto
    quality_format: standard
    strandedness: auto
    transcriptome_fasta: /swiftcache/ngs/oncology/analysis/todd_creasy/working/reference/GRCh38.primary_assembly.genome_GPC3CART.fa
    transcriptome_gtf: /swiftcache/ngs/oncology/analysis/todd_creasy/working/reference/gencode.v29.annotation_GPC3CART.gtf
  analysis: RNA-seq
  description: 12_LI6612_Baseline_8770
  files:
  - /swiftcache/ngs/oncology/analysis/todd_creasy/working/fastq/SM-745/SM-745-12_R1_001.fastq.gz
  - /swiftcache/ngs/oncology/analysis/todd_creasy/working/fastq/SM-745/SM-745-12_R2_001.fastq.gz
  genome_build: hg38
  metadata:
    panel: Baseline
fc_name: seqc
upload:
  dir: ../output/final

I recommend this is about Slurm rather than bcbio. I would personally do this as a Slurm job array (just my opinion). This is better to resolve via chat in my personal opinion. The CLI in the slurm script is key. I can see the Lmod via use and load probably via EasyBuild.  What I need to do is read about bcbio and see how thats put together.

You can see the bcbio docs for slurm here: https://bcbio-nextgen.readthedocs.io/en/latest/contents/parallel.html. The setup seems right according to the doc which is why I'm thinking there's an issue with the yaml

Okay, I've got it. One solution is an array job.  This will make it directly run in parallel rather than sequentially. Background this is EasyBuild on OpenStack (Swift) ...The GitHub request is here: https://github.com/bcbio/bcbio-nextgen/issues/3719

@M__ - I reached out to a bcbio expert and it was just a bad sbatch script. You can see his response here: https://github.com/bcbio/bcbio-nextgen/issues/3719

I understand (see Slurm array job solution below), what I can definitely help with is the automated generation of the seqc.yaml file. This avoids the repetition (and risk of error) in copy and pasting. The following one-liner exactly reproduces the seqc.yaml file and needs the chunk.yaml below.
for ((x=10;x<=13;x++)); do while IFS='' read -r a; do echo "${a//SM-745-10/SM-745-${x}}"; done < chunk.yaml >> seqc.yaml; done

This could be generalised to 3 command lines in bash,
y='SM-745-10'
z='M-745'
for ((x=10;x<=13;x++)); \
do while IFS='' read -r a; do echo "${a//${y}/${z}${x}}";\
done < chunk.yaml >> seqc.yaml; done

Here everything from for is a bash one-liner. Obviously ((x=10,x<=13;x++)) would be ((x=1,x<=20;x++)) for 20 fastq forward and reverse files. This will ensure the separate files are correctly named in the yaml file, in this example, viz.
# first pair
SM-745-10_R1_001.fastq.gz
SM-745-10_R2_001.fastq.gz
# second pair
SM-745-11_R1_001.fastq.gz
SM-745-11_R2_001.fastq.gz
# third pair
SM-745-12_R1_001.fastq.gz
SM-745-12_R2_001.fastq.gz

The chunk.yaml is,
details:
- algorithm:
    aligner: star
    disambiguate: mm10
    expression_caller:
    - salmon
    - kallisto
    quality_format: standard
    strandedness: auto
    transcriptome_fasta: /swiftcache/ngs/oncology/analysis/todd_creasy/working/reference/GRCh38.primary_assembly.genome_GPC3CART.fa
    transcriptome_gtf: /swiftcache/ngs/oncology/analysis/todd_creasy/working/reference/gencode.v29.annotation_GPC3CART.gtf
  analysis: RNA-seq
  description: 10_LI6612_Baseline_8743
  files:
  - /swiftcache/ngs/oncology/analysis/todd_creasy/working/fastq/SM-745/SM-745-10_R1_001.fastq.gz
  - /swiftcache/ngs/oncology/analysis/todd_creasy/working/fastq/SM-745/SM-745-10_R2_001.fastq.gz
  genome_build: hg38
  metadata:
    panel: Baseline
-

bcbio
To summarise the response from the Git support team is there's conflict occurring  ...

The issue is here for the slurm script ...

#SBATCH -p core -n 32

I assume this is -N 1 ,i.e. >= 32 core node, if it was 16 cores per node then -N 2.
However, bcbio is taking care of this and the restriction 32 is conflicting with 72 from bcbio, i.e. it creates its own Slurm controller script.

in bcbio yaml

    default:
    memory: 30G
    cores: 72

2a. Firstly 72 is specified in Slurm's CLI here ... -n 72
2b. This isn't needed because,
#SBATCH --mem-per-cpu=40G

--mem-per-cpu  This is a very recent feature BTW.

Slurm array job solution
What I considered was using the chunk.yaml above with an array job - this solution isn't needed because bcbio is dynamically reading and parallelising seqc.yaml (maybe deploying a wrapper for sinfo, squeue). In the solution below what I'm doing is placing each fastq pair directly into the Slurm queue. This would then parallelise everything directly. Again whilst it would work, it's unnecessary.
Thus for all 20 fastq files assuming SM-745-1_R1_001.fastq.gz ... SM-745-20_R1_001.fastq.gz.
#!/bin/bash

#SBATCH --job-name=bcbiopipeline_GPC3_PDx
#SBATCH -t 1-00:00:00 # 1- is one day. 00:00:00 is a wall clock. Does "1-00:00" work?
#SBATCH -c 1
#SBATCH --array=1-20 # this is now in $SLURM_ARRAY_TASK_ID
#SBATCH --mem-per-cpu=40G
#SBATCH --error=job.%J.err
#SBATCH --output=job.%J.out
#SBATCH --export=ALL

echo "Load bcbio..."
module use /projects/ngs_prodifx/local/software/modules/
module load bcbio-nextgen/1.2.9

# SM-745-10 and SM-745 could be assigned as variables y and z as above

while IFS='' read -r a; \
do echo "${a//SM-745-10/SM-745-$SLURM_ARRAY_TASK_ID"; \
done < chunk.yaml > seqc.yaml; done
bcbio_nextgen.py seqc.yaml -n 72 -t ipython -s slurm -q core -r t=0-72:00 -r conmem=20 --timeout 3000 --tag "rna"

Note the -n 72 parameter, i.e. number of cores deployed, will likely be reduced under this solution.
Alternative The alternative is to issue a slight modification of the bash one-liner described above to produce 3 (or 20) separate yaml files via bash while referring the different fastq paired files and use array job to cycle through them and delete each in turn when done.


text:
What happened to the PhyLoTA browser?
A paper I'm reading from 2021[1] mentions doing some analysis on data from the PhyLoTA database, with this reference citation:

M. J. Sanderson et. al. 2008 The PhyLoTa Browser: processing GenBank for molecular phylogenetics research, Syst. Biol. 57(3) 335-346.

This paper abstract claims to have an online interface at http://loco.biosci.arizona.edu/pb, but this link is not currently responsive. Does anyone know what happened to this project, and if it is currently accessible at some other location?

[1] S. M. Harrington et. al. 2021 Properties of Markov chain Monte Carlo performance across many empirial alignments MBE

Phylota was implemented in R and is present on GitHub here, described by Bennett et al (2018).
Here's the code...
library(phylotaR)
wd <- '[FILEPATH TO WORKING DIRECTORY]'
ncbi_dr <- '[FILEPATH TO COMPILED BLAST+ TOOLS]'
txid <- 9443  # primates ID
setup(wd = wd, txid = txid, ncbi_dr = ncbi_dr)
run(wd = wd)

Sanderson et al  was 15 years ago, albeit Sys. Biol. is a very reputable journal in phylogenetics. I've deployed solutions published there with outstanding results.


text:
Convert Symmetrix Matrix into Lower Triangular Matrix
I'm using Scipy's pdist function to calculate a distance matrix which is working great but it gives me a symmetric matrix and I want to convert it into a lower-left triangular matrix.
This is what it prints out for me:
  a   b   c   d
a 0 ### ### ###
b ### 0 ### ###
c ### ### 0 ###
d ### ### ### 0

This is what I want it to be:
###
### ###
### ### ###

I know there's a function in R to do this (upper.tri) but it just replaces everything to the right of the 0 with more 0s. Is there a better/easier way to get it like this? I'm opening the file in excel and removing everything to the right with the 0 manually right now.

From reading the scipy documentation for the pdist function, it seems to suggest it should already only return upper-triangular entries rather than a full matrix. Are you using pdist directly or through some other package?

This is the command I use to calculate the distance: pd.DataFrame(squareform(pdist(kmer_matrix, 'euclidean'))) I convert it into a pandas dataframe so it's easier to output into a file and I use squareform because I read somewhere that it helps space out the final output.

I just tried running it without squareform and it outputs all the distances in a descending list :/

NumPy has a convenient built-in function for accessing the lower-triangular or upper-trianglar entries of an array.
To get the lower-triangular entries, you can use numpy.tril_indices
to get the indices of these entries, and then use these indices to extract the entries.
import numpy as np

ar = np.array(
    [[0, 1, 2],
    [3, 4, 5],
    [6, 7, 8]]
)
idx = np.tril_indices(3, -1)
ar[idx]
# outputs array([3, 6, 7])

This doesn't automatically produce the spacing you want, but formatting the final display is likely context-dependent.

This could work for now thank you!!!

There seems to be a scipy method that retains the lower triangle but with the 0s. See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.tril.html#scipy.sparse.tril

Thanks! This is pretty helpful too


text:
BioPython SeqIO fails with Invalid Alphabet found
I am trying to store a GenBank file that could be read by Geneious later.
I found an example how to create the GenBank file here
But now I run into Alphabet issues.
I tried all protein classes from the documentation, but all fail with the same error:
TypeError: Invalid alphabet found, <class 'Bio.Alphabet.ProteinAlphabet'>.
I dont know what else to do here? Using biopython==1.73
    record = SeqRecord.SeqRecord(
        Seq.Seq("MIRQALAVAALLLAGTAQADGLIDN", alphabet=Alphabet.ProteinAlphabet),
        id="FooBar",
    )

    with NamedTemporaryFile(mode="w") as genbank:
        SeqIO.write(record, genbank, "genbank")
        pass

If possible, please show a minimum reproducible chunk of code. The code shown here does not include any library imports, so will not work as-is.

Quickly ...
from Bio.Alphabet import generic_protein
from Bio.Seq import Seq

then ..
myprot = Seq('EXAMPLEPROT', alphabet=generic_protein)

Thanks, I tried without setting the alphabet but Seq needs that parameter. Let me try another way of importing the Alphabet

Using generic_alphabet tells me this: ValueError: Need a Nucleotide or Protein alphabet

I was wrong, generic_protein does the trick, not generic_alphabet

Thanks @ElDude done .. many thanks.


text:
How can I resolve `bind_rows()` error?
So I wanted to download massive amounts of COX1 (or COI) sequences from any available database. I found this script in https://bioinformatics.stackexchange.com/a/13187/16657 this comment that achieves just that using the bold R package. I changed the script to fit my data preference:
# load packages
library(tidyverse)
library(rentrez)
library(bold)    # API interface to BOLD
library(taxize)  # for NCBI taxonomy lookup
library(seqinr)  # for FASTA output

set_entrez_key()

# get class-level taxa within "Mollusca" from NCBI taxonomy
taxa <- downstream("Mollusca", db = "ncbi", downto = "order")
# check if taxa present in BOLD
checks <- bold_tax_name(taxa$Mollusca$childtaxa_name)
taxa_bold <- checks[!is.na(checks$taxon),]$taxon

# Download sequences from BOLD for each class-level taxon
sequences <- map(taxa_bold, bold_seq, marker = 'COI-5P') %>%
  flatten() %>%
  bind_rows()

# Convert the list of sequences to a data frame
sequences_df <- do.call(rbind, sequences)

# Write sequences to a file
write.fasta(
  sequences = as.list(sequences_df$sequence),
  names = as.list(sequences_df$id),
  nbchar = 80,
  file.out = 'coi5p.fasta'
)

In the aforementioned thread, this script worked but now that I ran it I got the following error:
Error in `bind_rows()`:
! Argument 1 must be a data frame or a named atomic vector.
Run `rlang::last_trace()` to see where the error occurred.

I am not sure what to do for this. I commented on the original thread but since it is kind of old I decided to post about it here.
Note: the set_entrez_key() requires a NCBI API key which you can get by registering on NCBI.
Edit #1: The sequences object before the flatten() command is a large list containing 53 elements. After the command, the object becomes a large list with 257 elements.

Can you please show the object as it appears before and after the 'flatten' command? It's unusual for me to see that command in R code, and this is just before the bind_rows command that you're concerned about, so it will probably be important for working out what's going wrong.

@gringer I edited the question with the details about the object.

Please show the head of the object, or the first few lines of output. Your text description is not detailed enough for me to understand what's going on.

I managed to find the solution. After @gringer mentioned that they had not seen the flatten command I decided to dig deeper. Turns out the command has been superseded  in purrr 1.0.0. So I will be posting the new script for anyone who seeks an updated version of the one in the old thread.
# load packages
library(tidyverse)
library(rentrez)
library(bold)    # API interface to BOLD
library(taxize)  # for NCBI taxonomy lookup
library(seqinr)  # for FASTA output
library(rBLAST)

set_entrez_key("d00fee673c4499df4c2479ba5fd71850b308") #my ncbi api key, this is essential for pulling the taxonomy through ncbi

# get class-level taxa within "Mollusca" from NCBI taxonomy
taxa <- downstream("Mollusca", db = "ncbi", downto = "class") #setting the taxonomic level of the search. it will pull sequences for all mollusca but seperately for each class. if the search is large change to a lower taxonomic level

# check if taxa present in BOLD
checks <- bold_tax_name(taxa$Mollusca$childtaxa_name)
taxa_bold <- checks[!is.na(checks$taxon),]$taxon

# Download sequences from BOLD for each class-level taxon
sequences <- map(taxa_bold, bold_seq, marker = 'COI-5P')
# Convert the list of sequences to a data frame
sequences_df <- do.call(rbind, sequences)

# Filter sequences based on the "COI-5P" marker. Despite the markes parameter, some files may have integrated multiple loci. The other sequences are removed from the R object
coi5p_sequences <- sequences_df %>%
  filter(marker == "COI-5P")

# Write sequences to a file
write.fasta(
  sequences = as.list(coi5p_sequences$sequence),
  names = as.list(coi5p_sequences$id),
  nbchar = 60,
  file.out = 'coi5p_Mollusca.fasta'
)


text:
How to make a POST "vep/:species/region" request that is able to accomodate multiple alternate alleles for a single variant?
There is an R code example on Ensembl's website on how to execute a POST "vep/:species/region" request for variants that only have one alternate allele. The part that holds the variant data in the POST request is structured as follows:
..., body = '{ "variants" : ["21  26960070  rs116645811 G A . . .", "21  26965148  rs1135638 G A . . ." ] }'...

Where, for the first variant, 21 is the chromosome number, 26960070 is the start point coordinate on the chromosome, rs116645811 is the variant identifier, "G" is the reference allele and "A" is the alternate allele.
Alternatively, you can substitute the variant identifiers - rs116645811 and rs1135638, in the example above - for the endpoint coordinate of the variant. I use this option myself and can personally confirm that it works.
As mentioned before, the example above shows a POST request of variants that have only single alternate alleles.
Does anybody, on the other hand, know how (one should adjust the structure in the example above) to execute a POST request that accommodates multiple alternate alleles for a single variant?
If you explain your answer with an example like the one above, please use an example with the endpoint coordinate (and not the variant ID).
Thanks!

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Thanks for your question. The question has received a "closed" request. The site convention is for the OP to present an example for others - to  help them answer - here you as OP are  requesting the person answering to present the answer. This is can lead to misunderstandings about what is required. There are some very skilled R coders on the site who could probably answer this.


text:
Problem with edit and align viral whole genomes
I have a database of whole genome of African swine fever virus (about 190kb/genome) downloaded from GenBank. But these sequences do not start at the same place on the genome. When I align the sequence to build a phylogenetic tree. The sequences at the 3' and 5' ends have many gaps and do not match. The length of the entire data set is over 220kb due to many gaps.
Kindly suggest me a solution to the problem.

Welcome! Can you post an image of your alignment?


text:
Generating simulated gold standard NGS dataset/ Any publicly available smal gold standard NGS dataset
I am trying to check my variant generation pipeline for humans. The pipeline is for NGS. I am searching for any dataset(.fastq+.bam+.vcf) for any specific sample. Now, GIAB has a nice dataset for that. However, the problem is that it's not a tiny dataset. And it's not feasible to test my pipeline with this data every time. I tried to generate a synthetic dataset with some read simulators(NEAT, MASON2). However, their aligned read files(.bam) files seem not very good. You can check my Issue for NEAT in link.
So, my question is, Is there any publicly available tiny dataset(.fastq,.bam, .vcf) with which I can validate my pipeline fast? Or, is there any simulator with which I can simulate .bam, .fastq files from a given input .vcf? I am using fastqc,bwa,samtools,GATK4 mainly and GRCH38 reference in my pipeline.

I would suggest subsampling an existing dataset such as the GIAB resource that you've already indicated.
Use a utility such as seqtk to subsample the appropriate number of reads from the dataset, and use that.
Example with seqtk to get 1 million reads each from paired fastqs while preserving pairing:
  # note -s100 random seed which should be used for both to keep pairs
  seqtk sample -s100 giab_read1.fq 1000000 > giab_sub_r1.fq
  seqtk sample -s100 giab_read2.fq 1000000 > giab_sub_r2.fq

Ok, that's actually a good suggestion. However, how to subsample gold standard fastq, .bam, and .vcf? So, that I know my pipeline steps/changes are okay in all three steps?

@ShafayetRahat just subsample the fastq, and then run the small fastq files through your pipeline to get bam abd vcf. I use a similar approach, and it only takes about 5 minutes to run through the pipeline with these small files.

@ShafayetRahat terdon is right, unless you're doing something odd mapping should be unaffected by depth. And you will not be able to subsample a VCF to only represent certain reads, I am not sure how you would do that. (though in principle you can do it with a BAM file, it would just be a lot of work compared to re-mapping).


text:
How to identify TF binding motifs?
I have DNA sequences from a single gecko species for a set of genes (SULF1, SOX9, SAL1), I wanted to compare the putative promoter sequences for each gene, and compare them to the promoter regions of other related species. Would something like HOMER work for say finding TF binding motifs in around 20 different sequences for 20 different species, or would I need a larger dataset?


text:
Current 'most comprehensive' whole genome alignment resource for human?
I found the human genome alignment sets here:
https://www.ensembl.org/info/genome/compara/multiple_genome_alignments.html
but they doesn't include, for example, ancient human genomes, e.g. those here:
https://reich.hms.harvard.edu/allen-ancient-dna-resource-aadr-downloadable-genotypes-present-day-and-ancient-dna-data
or, for example, any of the other 'pan genome' whole genome alignments here:
https://humanpangenome.org/
However, those latter two resources don't include whole genome alignments...
About 10 years ago I remember reading a paper that used 'lots' of genomes to uncover novel features using different patterns of conservation, demonstrating the power of the approach.
Have I missed a new, more up to date, whole genome alignment resource?
Many thanks,

This seems like it might be an XY problem. Can you please provide more context to your question and explain what you want to do with this resource? You're asking for a fairly large dataset, and there's probably a better way to do what you want without downloading all that data.

The obvious thing missing from your list is the telomere-to-telomere genome assemblies that have been coming out in the last couple years.
Assembly link here.
There are now a few others as well that are similar.
It's been demonstrated that this T2T assembly improves mapping rates, in some cases pretty dramatically.
I am not sure entirely what you want from whole genome alignments, but UCSC if not others has incorporated the T2T assembly. Here is a related commentary (which I haven't read in any detail).