text:
Can someone help me estimating the runtime of the pipeline applied by the vertebrate genome project?
The vertebrate genome project (VGP)
has a lot of interesting publications such as this one.
The rough pipeline is outlined below:

Here the pipeline in more detail:

While the paper describes all the steps of the iterative assembly pipeline, I cannot find information on the runtime of their assembly-pipeline.
I was excited to check out their github page to find out more in the documentation. Unfortunately did I not find any information on the runtime in the https://github.com/VGP/vgp-assembly repo.
I'have written some e-mails asking collaborators of VGP but I have still not gotten any answers.
Did any of you apply this assembly-pipeline?
Could anyone tell me the runtime of the assembly-pipeline?
I hope this question is in the appropriate stackexchange, if not please point me to what you think fits best. Thanks for your suggestions, please be gentle I'm new to all of this.
For more info on the VGP check out this link

Just a comment (please don't turn this into an answer): like much work in genomics, the VGP paper is already outdated when it is published. These days, 10x is discontinued and few are using CLR. When budget allows, go for HiFi. HiFi assembly is much faster (half a day for mammals) and is of much higher quality.

Could you please specify - like posting a publication and a brand of a platform you'd recommend? Thanks in advance

From what I understand the VPG pipeline includes HiFi on the PacBio plattform anyways - so I'm not sure if you mean something else and newer...

Sorry for the late response. I haven't seen their v2.0 pipeline. Its paper, many repos and your original figure are about obsolete data types. For PacBio HiFi, you mostly need an assembler only. No need to apply the VGP pipeline.

@172818 Looking at papers I've seen only publications using a 'hybrid' approach... You say if I was to use PacBio HiFi I need an 'assembler only'. What workflow would you suggest? If I was to produce a high quality assembly, should I still use several platforms e.g. including illumina for polishing?

Note to OP that your tag of @user172818 was not successful in your last comment. That's probably why they didn't respond.
Because no one has gone there so far, I will simply link to HiFiAsm. See linked Nature Methods paper. This tool accepts HiFi reads, which are now standard, and gives you a GFA that you can quite easily convert to FA.
Polishing with Illumina would actively harm the assembly, as it would remove real heterozygosity or add errors where Illumina reads map ambiguously.
Run "out of the box"t with default options on HiFi it will probably give you much better results than VGP, haplotype-resolved.
If you also add Hi-C PE reads with no additional options, it will give you something much much better than VGP.
It's great. See the docs for practical info on running it.
There is no reason to follow VGP workflows, as they are very very outdated, as pointed out in comments.


text:
How to Read SCF file in Python?
Is there any way that I can read SCF file in python like in R using sangerseqR, I have tried with Biopython, it seems it does not support this format.

You could try converting it to ABI format as asked in this question and then import that

Alternatively, bioconvert supports conversion of SCF to FASTA and FASTQ formats.

Conversion to SCF to FASTA and FASTQ formats will not help me since I need the chromatography data not just the seq

It was done using sangerseqR with rpy2 in python
    import rpy2.robjects as r
    from rpy2.robjects.packages import importr
    utils = importr('utils')
    utils.install_packages('sangerseqR', repos="https://git.bioconductor.org/packages/sangerseqR")
    utils.chooseBioCmirror(ind=1) # select the first mirror in the list
    utils.install_packages('sangerseqR')

    sangerseqR = importr('sangerseqR')
    Trace = sangerseqR.readsangerseq( file="1I1_F_P1815443_047.scf") #scf sequence file is now loaded

sangerseqR Documentation(pdf)

source: Stack Exchange

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

sort: votes


text:
Different results in differential expression analysis of microarray data
I am performing differential gene expression analysis to microarray data for type 2 diabetes donors and nondiabetic donors. When I run the code I get some different results in each time (about 50 or so different genes in each run). Is this normal? I am submitting a project and it could be very problematic if my results change in each run, since I will perform downstream analysis for the DE genes and use web tools such as DAVID for further analysis.
Is there a way to keep the results from changing in each run?
g <- getGEO("GSE20966",GSEMatrix = TRUE)
g<- g[[1]]
exprs(g)
exprs(g) <- log2(exprs(g))

#Diffrential gene exp

library(limma)
design <- model.matrix(~0+sampleInfo$characteristics_ch1.1)

colnames(design) <- c("NonDiabetic","T2D")

summary(exprs(g))

exp <- na.omit(exprs(g))
cutoff <- median(exp)

## TRUE or FALSE for whether each gene is "expressed" in each sample
is_expressed <- exprs(g) > cutoff

## Identify genes expressed in more than 2 samples

keep <- rowSums(is_expressed) > 2

## check how many genes are removed / retained.
table(keep)
keep<- na.omit(keep)

## subset to just those expressed genes
g<- g[keep,]
fit <- lmFit(exp, design)

head(fit$coefficients)
contrasts <- makeContrasts(T2D - NonDiabetic, levels=design)
fit2 <- contrasts.fit(fit, contrasts)
fit2 <- eBayes(fit2)
topTable(fit2)
topTable(fit2, coef=1)
decideTests(fit2)
table(decideTests(fit2))
aw <- arrayWeights(exp,design)
aw
fit <- lmFit(exp, design,
            weights = aw)
contrasts <- makeContrasts(T2D - NonDiabetic, levels=design)
fit2 <- contrasts.fit(fit, contrasts)
fit2 <- eBayes(fit2)
g <- na.omit(g)
anno <- fData(g)
anno
anno <- select(anno,`Gene Symbol`,ENTREZ_GENE_ID)
fit2$genes <- anno
topTable(fit2)
full_results <- topTable(fit2, number=Inf)
full_results <- tibble::rownames_to_column(full_results,"ID")
library(ggplot2)
dev.off()
ggplot(full_results,aes(x = logFC, y=B)) + geom_point()
p_cutoff <- 0.05
fc_cutoff <- 1

full_results %>%
  mutate(Significant = adj.P.Val < p_cutoff, abs(logFC) > fc_cutoff )
 %>%
  mutate(Rank = 1:n(), Label = ifelse(Rank < topN, anno$`Gene
 Symbol`,"")) %>%
  ggplot(aes(x = logFC, y = B, col=Significant,label=Label)) +
geom_point() + geom_text_repel(col="black")

 p_cutoff <- 0.05
 fc_cutoff <- 1

filter(full_results, adj.P.Val < 0.05, abs(logFC) > 1)
 library(pheatmap)
topN <- 20

top<-150
gene_names <- mutate(full_results, Rank = 1:n()) %>%
  filter(Rank < top) %>%
  pull(Gene.Symbol)

DE results do not change from one run to another. The code given in your question will give identical results each time you run it.

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


text:
WGCNA hub genes error
I am trying to find hub genes for modules made using WGCNA:
hub = chooseTopHubInEachModule(datExpr, moduleColors)

But I am getting error:
Error in hubs[m] <- colnames(adj)[hub] : replacement has length zero

Thank you for the help1

Cross-posted on biostars: https://www.biostars.org/p/9486378/ Please don't shotgun post on multiple fora. See this discussion: https://bioinformatics.meta.stackexchange.com/questions/265/about-posting-queries-in-stack-bioinformatics

Can you please share the code, how you fixed this problem? I am facing the same problem and my data has colnames.
Any help will be appreciated.
Thanks

The most likely culprit here is that your expression data (the datExpr input to chooseTopHubInEachModule) has no colnames. Check that your data do have appropriate column names. Also check that the module color vector you supply has at least actual module (type table(moduleColors) and check that you have at least one module apart from "grey" or 0).

thanks as it worked. I have another question for extracting genes within modules of interest. I used unsigned network and interested in two modules for a trait. Earlier I made scatter plot and extracted the genes in all the modules. But now I only want to extract genes in these two modules for further analysis. I am not able to find a tutorial with more than one module. Can you please suggest something? Thank you1

Well, there are a few sections of the tutorial that may help. You can export the entire network analysis results (which you should do anyway) as illustrated in section 3d and then restrict the table (or just the gene identifiers) to the two modules you want. Alternatively, section 4a shows how to export gene identifiers in selected modules.


text:
Counting the number of gene isoforms in a GFF3, is this method correct?
Recently I've been tasked to count the numbers of gene isoforms for each locus in a .gff3 file. I'm still doing my first steps in biology and bioinformatics, so I still struggle with some concepts and methodologies.
According to my understanding a gene isoform is a mRNA located in the same locus but with different start positions, so I did the following procedures:

Selected the entries in the .gff3 that are classified as mRNA.
In my .gff3 there is the additional information of the locus of the mRNA, so I checked for each locus how many unique values there were for the Start values and subtracted one.
Assigned each value from step 2. to their respective locus, resulting in the desired number of isoforms for each locus.

My question is is this method correct?
Why do I believe it could be wrong or inaccurate:
I don't have a complete solid understanding of a gene isoform, and researching some more I've seen that not only the start point of the mRNA is relevant but other things such as the protein coding DNA sequence (which I believe could be identified in the file under the CDS tag) and untranslated regions (which I don't know if I can identify in my file), and in the method I described above I've only taken in consideration the Start position, what appears to be incomplete.
Additionally if the method is incorrect how can I take in consideration the additional circumstances, and if it's correct why  didn't those circumstances influence the count so much?

In my .gff3 there is the additional information of the locus of the mRNA, so I checked for each locus how many unique values there were for the Start values and subtracted one;

Different isoforms of mRNA can share the same start coordinate. For example, they could share the same start coordinate but one skips an exon. MSpaint illustration below to help. Generally for GFF3 I would simply count e.g. how many mRNA have Parent=GENE0001234 to get the number of mRNA for GENE0001234, and you can expand that for all instances of GENENNNNNNNN

See also for example Wikipedia for Alternative splicing https://en.wikipedia.org/wiki/Alternative_splicing

Gene isoforms are transcripts derived from the same gene that vary in their exons (blocks of the gene that make up the final transcript). This can alter the start position, end position, and/or internal sequence of the transcript. GTF/GFF files account for this by having multiple transcripts per gene (if applicable) that capture the various isoforms of the gene.
Here's a little example in R to count the number of transcripts for a few human genes found in an ENSEMBL reference GFF3 file.
library("GenomicFeatures")

gff3 <- makeTxDbFromGFF("http://ftp.ensembl.org/pub/release-104/gff3/homo_sapiens/Homo_sapiens.GRCh38.104.chr.gff3.gz")

> stack(vapply(transcriptsBy(gff3, "gene")[1:10], length, 1L))
   values             ind
1       5 ENSG00000000003
2       2 ENSG00000000005
3      16 ENSG00000000419
4       5 ENSG00000000457
5       9 ENSG00000000460
6       7 ENSG00000000938
7       6 ENSG00000000971
8       3 ENSG00000001036
9      14 ENSG00000001084
10      2 ENSG00000001167


text:
GROMACS RMSF Result Analysis Problem
This is RMSF for protein from GROMACS after 100ns simulation-

I have few questions-

Residue numbers are starting from 1650 to 1850, has it ignore LIG
and HOH as residues?
There are three lines, so are there three chains and their residue number is same?
Are the straight lines connecting these chains or this is kind of loop something?
Since I am seeing 2nm or 20A fluctuation, is the protein quite disturbed?

The residue numbering in a PDB file mainly depend on the person who deposited it. There are usually reasons for how things are numbered, but ultimately you should look at PDB file, and check the details on how the structure was obtained. To give you an example, during my PhD I worked on a protein that had some missing residues at the N-terminus. That was because the N-term was very unstable and they didn't manage to crystallise it properly.
You need to check the PDB, but it is possible that the 3 chains have the same residue numbering (A1-200, B1-200). This is common in structures where you have multiple identical sub-units (e.g. dimers).
The lines are merely an artefact. Those lines are the attempt of the engine to draw a single line for all the chains. You go from residue 1850 of chain 1 to residue 1650 of chain 2, and so on. Some tools handle this behaviour, others don't. I would suggest you to parse the PDB yourself and plot (e.g. ggplot) each chain as a layer of the same figures .
It is indeed, but only you can tell if this is good or bad. It really depends on how you calculate the RMSF. Was it on the backbone or on the side-chains? How is your dynamics? High temperature? Does you protein partially unfold? are the chains separating?

Thanks for your response and clarification to understand the result.


text:
BWA-MEM2 alignment-Snakemake
I have started using snakemake 6.5.2 to align fastq files with reference file. I have pasted the error below in this question. How to allocate memory in the snakefile and read the header from samfile, '-'. This is the snakefile (wrapper for running alignment):
 rule bwa_mem2_mem:
        input:
            reads=["/scicore/home/cichon/GROUP/test_workflow/samples/{sample}.1.fq", "/scicore/home/cichon/GROUP/test_workflow/samples/{sample}.2.fq"]
        output:
            "/scicore/home/cichon/GROUP/test_workflow/output/{sample}.bam"
        log:
            "/scicore/home/cichon/GROUP/test_workflow/output/{sample}.log"
        params:
            index="/scicore/home/cichon/GROUP/test_workflow/data/gch38.fa",
            extra=r"-R '@RG\tID:{sample}\tSM:{sample}'",
            sort="none",             # Can be 'none', 'samtools' or 'picard'.
            sort_order="coordinate", # Can be 'coordinate' (default) or 'queryname'.
            sort_extra=""

        threads: 8
        wrapper:
            "0.77.0/bio/bwa-mem2/mem"

I used this command to execute the above snakemake file:
snakemake --use-conda --snakefile Snakefile --cores -p /scicore/home/cichon/GROUP/test_workflow/output/DE98NGSUKBD117612_1.bam /scicore/home/cichon/GROUP/test_workflow/output/DE71NGSUKBD117613_3.bam

Error from log file:
Looking to launch executable "/scicore/home/cichon/GROUP/test_workflow/.snakemake/conda/0ecf6ba401c70d7501c73bb1c733a060/bin/bwa-mem2.avx", simd = .avx
Launching executable "/scicore/home/cichon/GROUP/test_workflow/.snakemake/conda/0ecf6ba401c70d7501c73bb1c733a060/bin/bwa-mem2.avx"
-----------------------------
Executing in AVX mode!!
-----------------------------
* SA compression enabled with xfactor: 8
* Ref file: /scicore/home/cichon/GROUP/test_workflow/data/gch38.fa
* Entering FMI_search
* Index file found. Loading index from /scicore/home/cichon/GROUP/test_workflow/data/gch38.fa.bwt.2bit.64
* Reference seq len for bi-index = 6544178411
* sentinel-index: 2782712944
* Count:
0,  1
1,  1914743810
2,  3272089206
3,  4629434602
4,  6544178411

* Reading other elements of the index from files /scicore/home/cichon/GROUP/test_workflow/data/gch38.fa
* Index prefix: /scicore/home/cichon/GROUP/test_workflow/data/gch38.fa
* Read 0 ALT contigs
[fread] Cannot allocate memory
[main_samview] fail to read the header from "-".

Thanks
Priya

You could at least reference the biostars thread: https://www.biostars.org/p/9486284/
It is the same problem. Did you upgrade your snakemake and apply the suggestions? Voting to close this here to keep it focused in a single thread.

I’m voting to close this question because it is already being discussed at biostars extensively.

Please don't close the question. I am the one posted the question in Biostars as I didn't get a clear answer, I have posted it here again . Thanks

Did you follow the suggestions? If so please add comment and follow up rather than starting from scratch here.

Clearly, this seems a problem with memory allocation. You should tell us what you tried from the suggestions of Biostars so we can help you in your problem.


text:
dbSNP GRch37 Download
I am still new to this.
I am trying to download the GRch37 dbSNP data from the NCBI dbSNP Browser using the "Send to:" button. But whenever I select
"Send to:" > "File" > "Tab" > "Default order"
The coordinates of my positions will be GRch38 based. Does anyone knows how can I download the dbSNP data with GRch37 coordinates?
Thanks
Ken

I see no option to specify the reference version using the NCBI dbSNP browser, it is set to download the latest built information. I think this is not the standard tool to download dbSNP data across many sites, but rather to browse it and visualize it for specific variation.
If you are only interested in getting info for a smaller range of positions, I recommend you to use the UCSC browser, go to Table Browser, -> assembly GRch37, -> group variation, -> track all SNPs, -> table snp51 -> have a look at "describe table schema" to see what will you get, select your genomic range, and click on get output.
Alternatively and if you have a lot of space available in a server, simply go to the NCBI dbSNP FTP site, and download the latest GRch37 dataset (GCF_000001405.25) with all the entries. Make sure you read the README for further info about the files and data.


text:
Mapping unique GO term description given a specific GO ids
This question was also asked on Biostars
I have a list of GO ids and I want to find out a unique term description such that if I provide say 200 GO IDs I will give 200 specific GO terms. The code snippet I am using is given below. Although it kind of works, for one given  GO ID I get several terms like for my 200 GO ids I get 2577 terms whereas I am interested in a specific GO term. The output that I want is like this GO:0007254 = JNK cascade.
I will like to mention that I have also cross-posted the query in biostars and got no reply for the last 2 days and hence I am posting the query on Stack exchange bioinformatics. Any help will be useful.
library(biomaRT)
mart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL", dataset = "drerio_gene_ensembl", host = "asia.ensembl.org")
v1 <-c(target)   # The list of GO ids that I want to pass as a list.
go_list <- getBM(attributes=c("go_id", "name_1006", "namespace_1003"),
              filters = "go",
              values = c("v1"),
              mart=mart, uniqueRows = TRUE)

library(biomaRt)
mart <- biomaRt:::useMart(biomart = "ENSEMBL_MART_ENSEMBL", dataset = "drerio_gene_ensembl", host = "asia.ensembl.org")

v1 <- c("GO::0009055","GO:0016020","GO:0016491","GO:0006122","GO:0008121")

go_list <- biomaRt:::getBM(attributes=c("go_id", "name_1006", "namespace_1003"),
              filters = "go",
              values = v1,
              mart=mart, uniqueRows = TRUE)

so it is likely the "v1" that should not be in quotes
head(go_list)
        go_id                                                   name_1006
1 GO:0016020                                                    membrane
2 GO:0005739                                               mitochondrion
3 GO:0005743                                mitochondrial inner membrane
4 GO:0070469                                                 respirasome
5 GO:0006122 mitochondrial electron transport, ubiquinol to cytochrome c
6 GO:0001525                                                angiogenesis
      namespace_1003
1 cellular_component
2 cellular_component
3 cellular_component
4 cellular_component
5 biological_process
6 biological_process

when search filters option is used  searchFilters(ensembl, pattern = "GO") . We get the above filters (ensembl, pattern = "GO")
name                        description
12           with_go                      With GO ID(s)
13   with_goslim_goa              With GOSlim GOA ID(s)
58                go         GO ID(s) [e.g. GO:0000003]
59        goslim_goa GOSlim GOA ID(s) [e.g. GO:0000003]
183 go_evidence_code                   GO Evidence code. NO go_id filter exists.


text:
How to convert vcf to tbed/tfam and again convert it back to tbed/tfam to vcf using plink
When converting vcf to tbed/tfam format, I am getting error like
unrecognized flag<'--transpose'>

Code I used:
plink2 --vcf sample_vcf_file.vcf --recode --transpose --out outputfile

Correct command is
plink2 --vcf sample_vcf_file.vcf --recode transpose --out outputfile

vcftools also can convert vcf into tped/tfam

Yep, from plink manual: --transpose Deprecated. Use "--recode transpose"..


text:
Same region/region but different annotation
I have done the differential accessibility analysis then annotated the region with chipseeker.
my output for one of the many region and question regarding this
Peak log2FoldChange         padj seqnames    start      end width strand        annotation geneChr geneStart
105  Peak32946      -1.992474 9.197721e-19    chr11 14078815 14079558   744      * Distal Intergenic      11  13962689
5445 Peak32938      -2.473959 4.878995e-04    chr11 13961467 13961818   352      *          Promoter      11  13962689
      geneEnd geneLength geneStrand geneId      transcriptId distanceToTSS         ENSEMBL SYMBOL  GENENAME
105  14267884     305196          1  10418 ENST00000576479.3        116126 ENSG00000262655  SPON1 spondin 1
5445 14267884     305196          1  10418 ENST00000576479.3          -871 ENSG00000262655  SPON1 spondin 1

So the same region is marked as Distal Intergenic and Promoter. Next I have the same comparison RNA seq expression data. So now my objective is to find those genes+region which are expressed as well as deferentially accessible.
So for example here in this case in RNA seq data there is only one SPON1 gene so if I have to show it associated region how to put that,shall I consider promoter or the distal intergenic or consider which is more significant.

You have two different peaks at different positions with different annotations. This makes complete sense, you're just misreading it.

"you're just misreading it." this part I got it the gene got two different binding sites. So if i have to find if the accessibility correlates to the expression shall i consider both of them such as taking the mean ?

Does it make biological sense to average them? No, it does not. Which is biologically more likely to have an effect? The answer to that will tell you what to do.


text:
How can I calculate the distances between two specific residues of a protein from a PDB file?
Suppose, we have a protein (16PK).
We are considering a 5-residue segment/window.
If we assume ASN10 to be the ith residue, how can I calculate the following distances in Python?

from [i+2] to [i-2]
from [i+1] to [i-2]
from [i+2] to [i-1]

Next major question is, what would be the distances like if my segment/window is 7 or 9 residues long?

See this question for recommendations on a number of tools to calculate such distances: https://bioinformatics.stackexchange.com/questions/783/how-can-we-find-the-distance-between-all-residues-in-a-pdb-file

Your $[i,j,k]$ coordinates in the PDB file are coordinates. You can readily compute Euclidean distances between any pair of coordinates as described in the wiki page:
$d(p,q)=\sqrt{(p_1-q_1)^2 + (p_2-q_2)^2 + (p_3-q_3)^2}.
$
This resource suggests one such approach (code reproduced below in case that site stops existing):
import Bio.PDB
import numpy

def calc_residue_dist(residue_one, residue_two) :
    """Returns the C-alpha distance between two residues"""
    diff_vector  = residue_one["CA"].coord - residue_two["CA"].coord
    return numpy.sqrt(numpy.sum(diff_vector * diff_vector))


text:
How to get a consensus sequence from a nanopore fastq files?
I am new in bioinformatic field. I would like to know a way to generate a consensus sequence from nanopore fastq files (fastq files demultiplexed).
I usually generate a consensus sequence with "MEDAKA" of artic protocol, but it is used to generate a SARS CoV 2 sequence only.
I would like to know a differente way to do it with another program, R, or python.
Thank you!

When you say "consensus sequence", do you mean a genome sequence (genome assembly)? Or do you want to correct the reads? Medaka is normally used to correct existing assemblies. It looks like ARTIC is mostly amplicon-based, though I could be wrong about that. If what you want is to take fastq files and assemble a new genome sequence without a reference sequence, you probably want a genome assembler. Are your reads DNA, cDNA, direct RNA sequencing, something else? If RNA, maybe try rnaSPAdes: https://cab.spbu.ru/software/rnaspades/. If DNA, you could try raven: https://github.com/lbcb-sci/raven

I think you are confused about the process. Medaka is not covid specific. Use Medaka if you've used it in the past or are familiar with it.


text:
How to get a FASTA format file with the DNA sequences of all annotated genes?
I am analysing Pyrococcus Furiosus DNA sequencing data by considering data published here in NCBI. When I click on "Send to">"Gene Features">"FASTA format" I download a file that has the sequences of genes of this organism but I realized that this file has some sequences that are doubled (that is the sequence is present twice with different identifiers)... is there a way to get a file with all genes annotated and their respective DNA sequences without double sequences and with all sequences in the 5' to 3' direction without reverse complement in a well "ordered" way ? In NCBI (in the link I reported before) it indicates 2,128 genes so I would like a file with all these 2,128 genes annotated and their respective DNA sequences in FASTA format. Do you know if there is an other website or an other place in NCBI in which I can search to get this kind of file?
EDIT: if you follow the path "Send to">"Gene Features">"FASTA format", now it seems that there are just 2,128 genes so this problem is vanished in this case. Anyway, I would appreciate the answer for the future and for the reverse complement gene sequences that are still there.
Thank you very much.

From this SO post. You can use seqkit to remove duplicate sequences with the command below:
seqkit rmdup -s < sequences.txt > out.fa

The rmdup option removes duplicates, and the -s option calls duplicates on the basis of sequence, ignoring differences in headers.

thank you @Throckmorton I will try ! Even if my question was more focus on if it exists some website (in general) that has a FASTA format file with all annotated genes, without double DNA sequences or sequences that are reverse complement etc... so that if I remove the FASTA headers, I would obtain the complete genome DNA sequence

@Manuela does this species really have nothing but genes in its genome? And all of its genes are unique? In other words, are you sure that concatenating all genes will give you the genome? That seems very unlikely to me, but then I've spent most of my career working on eukaryotes.

@terdon thank you for the observation ! I assumed it was so... I am not a biologist :)

Ah, then @Manuela it is almost certainly the case the the genome has a lot of DNA that isn't part of any gene. I don't know archea at all, but if we take human as an example, genes are only around ~5% of the genome. If you want the genome, then look for the genome and not the genes.

ok @terdon thank you .

Sorry @terdon , but I have a doubt related to this: is it completely wrong to search a DNA sequence motif in a genome sequence  (because I could find a DNA sequence similar to my DNA sequence motif within the part that does not belong to gene and then I count it as a motif but it can not be) ?

@Manuela that depends entirely on what you are looking for, what species you are working on etc. If your motif has a function that is only relevant within a coding sequence, for example if you are looking for a DNA sequence which when translated produces a protein domain, then you want to focus on coding sequences only. If you want to see how common the motif is in general, then you can look at the whole genome. But there is no simple answer, it will always depend on what question you are asking.

ok thank you very much.


text:
Translate all reads in a .fastq into protein sequence from deep mutational scanning experiment
I'm working with paired-end NGS reads from an Illumina platform. The sequences I have are all of the same gene, but have one or more substitution mutations each.
Here is a rough workflow for generating the reads:

Expression of a library of sequence variants for a heterogeneous protein of interest in phage.
Selection of high-fitness variants according to a custom assay.
Extraction of the variant pool via amplification of the variable region of the gene of interest.
Preparation of Illumina sequencing library from the amplicon pool.
Deep sequencing of the amplicon pool.

I've already combined the paired-end reads with USearch and manually removed sequences with subpar quality scores, but I now need to translate all of the remaining DNA sequences into amino acid sequences. I'm uncertain of where to begin addressing this issue, and I was hoping someone could point me in the direction of an application or python package that could assist me.
My goal: Count the abundances of each sequence variant within the Illumina data to identify which of the protein sequences are the most fit.

Won't you need to map to a reference? You will need to know what frame to use to translate the reads. Unclear to me what the experimental design is here. Where do these reads come from? Is it an amplicon library? (If it is amplicon it is possible that you don't need to map anywhere.)

Reads are coming from a mixed population of phages which express multiple sequence variants of a protein of interest on their surface. The protein of interest is then put under selection. after selection the gene sequences of phages were collected and quantified to identify copy number of each variant in the population. This allows us to determine the fitness of each variant by seeing how many of that sequence there is.
I suppose I might need a reference map, but I am pretty inexperienced to know which program would be best.

"gene sequences of phages were collected" how? Presumably by sequencing, but what was the sequencing library preparation procedure? Was it whole genome shotgun e.g. nextera prep, locus-specific amplification (i.e. your protein of interest), barcode amplicon, something else? Presumably you limited sequencing to your protein of interest somehow, otherwise this would be a somewhat inefficient workflow.

Thank you for these questions! This is not my wet experiment, I am simply trying to use the data. But I can tell you the read libraries were generated with locus specific PCR.

Ok, one further question: Do all of your forward reads start at the same coordinate? E.g. are most reads mostly the same sequence? If paired-end, same question for reverse reads.

I already accepted your answer, but yes this is the case! Thank you very much!

I might edit your question slightly to include additional relevant info.

Based on back and forth in the comments, I think that you might be interested in a deep mutational scanning workflow, for example dms_tools2 or enrich2.
These are built specifically for the kind of experiment that you are interested in, if I understand correctly. See attached figure from dms_tools2:

Thank you! I think I will be using enrich2 for this data processing, it matches very well with my needs!


text:
no expression for hub genes identified in wgcna
I tried to identify intramodular hub genes in WGCNA analysis. I used:
hub = chooseTopHubInEachModule(datExpr, moduleColors)

But when I checked the hub genes for the modules, I found that the hub genes are expressed only in less than 5 out of 300 samples. I did not filter the gene expression data for any low expressed genes in certain number of samples. How does this make biological relevance if the hub genes are not expressed in most of the samples?
Thank you!

It could be that your modules group together genes (over)expressed in only certain small groups of samples that may be biologically relevant. You would need to check whether the samples in which the hubs (and presumably the modules) are up share some sort of a characteristic.
In general though, there are good reasons why authors of WGCNA recommend filtering out non- and low-expressed genes - correlations tend to not be a good co-expression measure when counts are very low, and technical effects can be much more pronounced. It is also recommended to spend a bit of time on preprocessing to make sure the data set does not have strong global drivers that are technical or uninteresting.

thank you! I applied some filters on RNA-seq data input, gsg = goodSamplesGenes(GE.adjusted,  minFraction = 1/2, minNSamples = 150, verbose = 3). I still have the genes which only have expression values for even one sample out of 300. I have a lot of zero's for no gene expression in the input file. Does WGCNA treat zero's as no or low gene expression? Should I replace zero's with NA so that it drop off the genes that have gene expression for less than say 100 samples? Thank you!


text:
How do I tell that a MinION run has finished based on the files produced?
I want to automate some scripts which should run after a MinION run completes. Is there a file which is produced by MinKNOW which will enable me to tell that the run has finished?

Use rampart and just stop it once you've happy with coverage

You can just look for the presence of a file called "final_summary_*.txt" in the run folder. From Matt Loose.


text:
Accessing expression data in an ExpressionSet
I used rma from the oligo package to normalize CEL files, returning an object of class ExpressionSet. How do I access the actual expression data?
> data = read.celfiles(names)
> data2 = oligo::rma(data)
class(data2)
[1] "ExpressionSet"
attr(,"package")
[1] "Biobase"
> data3 = as.matrix(data2)
Error in as.vector(data) :
  no method for coercing this S4 class to a vector

Have you looked at the oligo::rma documentation? Can you clarify why this doesn't make sense to you?

The class of your object which the rma function returns is an ExpressionSet which is a standardized data structure in Bioconductor/R from the Biobase package.
It stores both expression (a numeric matrix) data as well as column- and row annotation data.
You have to use the getter function exprs() to access the expression data slot, so exprs(data2) which is then of class matrix, link to documentation.


text:
How to get phylogenetic tree from multiple genes?
I constructed a phylogenetic tree using a gene (example - secA). I had to gather the same gene sequence for all the required species from public database-NCBI and then constructed the tree after multiple sequence alignment. I used the MEGA X software.
I now want to build a single phylogenetic tree using concatenated multiple genes, around 20.
How do I go about it.
Do I,

Concatenate all gene sequences for each species and do a multiple alignment or
perform an alignment for each gene and finally concatenate the alignments to get the tree?

Could this be done using MEGA X?


text:
java error when running CIRI-full.jar Merge module
I have an issue running a module trying to detect circRNA from RNAseq data. The error is a java error and I have not knowledge of java so I don't understand the issues and googling hasn't helped.
the issue arises when I try to run this module of the algorithm
https://ciri-cookbook.readthedocs.io/en/latest/CIRI-full.html#the-merge-module
as I try to run the module I get the following output and error
Merge module start

Loading Annotation...
Annotation Loaded
Loading CIRI_AS output...
CIRI_AS output Loaded
Loading CIRI output...
CIRI output Loaded
Loading CIRI_RO output...
CIRI_RO output Loaded
Combine AS and RO output
Combine completed, 71 reads are used
Exception in thread "main" java.lang.reflect.InvocationTargetException
        at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
        at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
        at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
        at java.lang.reflect.Method.invoke(Method.java:483)
        at org.eclipse.jdt.internal.jarinjarloader.JarRsrcLoader.main(JarRsrcLoader.java:58)
Caused by: java.lang.NullPointerException
        at Exon_gtf.get_exon(Exon_gtf.java:16)
        at Merge3.getmerge(Merge3.java:784)
        at CIRI_Full2.main(CIRI_Full2.java:551)
        ... 5 more

All prior steps preceding this one seem to have run properly and I was able to generate all the files to run the Merge module.
any help or explanation as to what the issue could be would be greatly appreciated. As I mentioned before I have no experience in java.
kind regards

Whatever the reason, it shouldn't crash with a null pointer exception.  You should bring that to the attention of the developers so they can try to debug.  Does it run successfully with their test dataset?  That might help narrow it down.

Hi Jesse, I did bring it up on github but no response so far. I have run the same pipeline of another dataset of mine and it worked. So the issue might be data I am feeding it. Thanks

Makes sense (something about the input data crashing the software).  It might be meaningful that the crash is within a GTF-handling function, too.  Are you giving it the optional GTF annotation file, and if so, does that same file work with the other dataset?  If you can make a small enough reproducible test case to post here that shows the crash, maybe someone can help narrow it down.

thanks for the feedback. I don't use the optional GTF. Comparing the inputs for that piece of code there is a clear difference in size so that means the difference is in there. I discussed with our collaborator who generated the data and they did process it differently compared to the previous one that did work. So the problem is not the code but the input and how the data was generated.


text:
Create BED file from list with gene symbols
Is there a way to generate a BED file of a list of gene symbols for a given genome?

Via Gencode and BEDOPS convert2bed:
$ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "gene"' - \
    | convert2bed -i gff - \
    > genes.bed

Replace the FTP path with that needed for your genome of interest. Feel free to edit your question if you need help with that, by specifying your given genome.
You can modify the awk statement to get exons, by replacing gene with exon.
BEDOPS: https://github.com/bedops/bedops
Once you have the genes in BED format, you can grep the file for them, using a text file of matches:
$ grep -Fwf genes_of_interest.txt genes.bed > genes_of_interest.bed

You could also get your genes from UCSC. Here's how:


text:
Alignment with inserts and keeping the indexing of ref seq intact
Parts of sequences are given below-
Reference sequence (pre-alignment):
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAACGAACTTTAAAATCTGTGT

Reference sequence (post-alignment) and below it is Sample sequence (post-alignment):
--------------------------------------------------------------------------------------------------attaaa---------ggtt------------------tataccttc---------ccaggtaacaaa-------------ccaacc-----aactttcgatctcttgtagatctgttctctaaacgaactttaaaatctgtgt
--------------------------------------------------------------------------------------------------------------------------------------------------------------taacaaa-------------ccaacc-----aactttcgatctcttgtagatctgttctctaaacgaactttaaaatctgtgt

I'm adding a simpler to interpret example-
Say the ref seq is aattaaatttgggggtttt and the sample seq is ttaaggggttaaatttgggggt--t. Then post-alignment, they will be like-
--aa----ttaaatttgggggtttt
ttaaggggttaaatttgggggt--t

Now, since my ref seq was
0    1    2    3    4    5    6    7    8    9    10   11   12   13   14   15   16   17   18
a    a    t    t    a    a    a    t    t    t    g    g    g    g    g    t    t    t    t

I want that post-alignment also, the indexing should be conserved-
           0    1                        2    3    4    5    6    7    8    9    10   11   12   13   14   15   16   17   18
 -    -    a    a    -    -    -    -    t    t    a    a    a    t    t    t    g    g    g    g    g    t    t    t    t
 t    t    a    a    g    g    g    g    t    t    a    a    a    t    t    t    g    g    g    g    g    t    -    -    t

I want to keep the indexing of the reference sequence conserved(i.e. the first base in ref seq post-alignment is a, second is t, third is t, etc.), like they do in standard softwares, and then I want to run some quick analysis on it, say to check for conservation of a (mono/di)-nucleotide at some positions. If anybody has some insight on how to do it the most efficient way(memory-wise and time-wise), then that'd be great. I use Python for my work.

May be missing something, but is there a reason you can't map the reference index to an alignment index using e.g. a list with the reference positions in the alignment? e.g. ref_posns = [3, 4, 6] if only the 4th, 5th, 7th positions are ungapped in reference; you can use this to access only those aligned positions. There is probably a more standard way to do it than that, though I don't know it off top of head. It seems rather complex to introduce an indexing scheme that doesn't have an index for many positions of the alignment. May help to know what "standard softwares" you refer to.


text:
Transfer annotations from one genome assembly to another
I've got an annotated draft genome assembly made using short and long read strategies. I've also done optical mapping to stitch some of the contigs together, with the goal being chromsome-length assemblies. How do I annotate the output .fasta from the optical mapping assembly using the .gff of the draft assembly that was used as input? I also have a .agp file that describes what sequence fragments ended up where in the output.

You should use a liftover approach. You could either do this based on exact changes, which might be a bit tedious given the presumably large number of join operations you've done (especially if you've e.g. broken contigs).
Therefore it might be more convenient to do it on the basis of homology.
liftoff is one simple tool for doing this with a gff file and your 2 fasta files. It seems to mostly work ok.


text:
How to see effect of point mutation in PyMOL?
I am trying to model how a specific point mutation would affect the protein structure.  So far, I've figured out how to create the point mutation in PyMOL using the mutagenesis wizard.  I'm not sure how to take the next step toward actually visualizing the consequences of the mutation though (i.e. would it disrupt a helix, etc).  Is there a way to do this in PyMOL, or would I need other software?  Thank you!

I would have said the question lacked detail, but the answers are good. Duet and equivalents are good. Pymol creates mutations, but can't analyse them.

PyMOL just find the rotamer that causes the least clashes, it does not repack other nearby residues, alter the backbone or adopt a conformer that is not part of the rotamer set it knows, this is a problem in some complicated cases core cases, but in in most cases the clashes reflect what the ∆∆G may be so I have to reluctantly agree here...

A common metric used to assess if a mutation destabilises the structure is the difference in Gibbs free energy. PyMOL does not give a ∆∆G (difference in Gibbs free energy, kcal/mol), but a very crude steric clash nor does it repack neighbouring sidechains out of the way, or alter the backbone, so unfortunately is not the way forward. However, if one must there is a way (see SIMBAI below).
Beyond destabilisation
Destabilisation predictors have a ~60% or lower success at predicting if something is pathogenic (ref). This is because the effect a mutation has on a protein is rather complicated as you can have destabilisation, but also aggregation, loss of regulation, mislocalisation, loss of binding to another protein, loss of binding to ligand, loss of activity (catalysis, inability to change conformation) etc etc. Then you can get curious subcases, such as loss of activity that result in the sequestration of a binding partner or homologues in a complex resulting in a worse phenotype than destabilisation or truncation.
Modelling
The behaviour of a protein can be assess by the energy of the system using forcefields, which are multi-term equations that describes the atomic interactions (attraction, repulsion etc.). The result is the Gibbs free energy mentioned above, which is a potential, so negative is good. Kinetic energy is release when the unfolded protein "rolls down" down the potential energy funnel —ProTherm dataset (an empirical dataset of ∆∆Gs) and the tool SDM give this value inverted as the negative can be confusing.
An uphill change in the Gibbs free energy (∆∆G > 0) generally correlates with a destabilisation of a protein —aggregation and solubility being the caveat in the "generally". The destabilisation a single state (static conformation) can be assessed. gmx energy in Gromacs, FoldX or Rosetta can be used for this.
For altered ability to switch conformation steered MD simulations are required, but generally the ∆∆G at different states can be used.
For the non-destabilisation options, https://www.phosphosite.org/ and http://elm.eu.org/ are very useful.
I have been working on a server that gives a ∆∆G value and note if there are any nearby post-translational modifications, gnomAD or changes in potential linear motifs etc.: Venus (unpublished). However, there are many ∆∆G calculation servers out there and are compared in many papers (example), but it should be said that the dataset used in these comparisons is always ultra-biased despite how much is tried. Missense by the Phyre2 group (Michael Sternberg) give detailed information of the neighbourhood changes the mutation causes.
In terms of extra info, Miscast gives Uniprot info.
SIMBA-I via PyMOL
In terms of PyMOL, the SIMBA-I algorithm can be used by getting the RSA from it —I would not recommend it as it may underestimate the severity of unusual cases. But if one must, in PyMOL one can paste the following into the console with the wild type variant:
# defined by you so change accordingly
from_residue = '' # single letter
to_residue = '' # single letter
resi =
chain = ''

## data
maxASA = {'A': 121.0, 'R': 265.0, 'N': 187.0, 'D': 187.0, 'C': 148.0, 'E': 214.0, 'Q': 214.0, 'G': 97.0, 'H': 216.0,\
          'I': 195.0, 'L': 191.0, 'K': 230.0, 'M': 203.0, 'F': 228.0, 'P': 154.0, 'S': 143.0, 'T': 163.0,\
              'W': 264.0, 'Y': 255.0, 'V': 165.0}
volumes = {'A': 88.6, 'R': 173.4, 'N': 114.1, 'D': 111.1, 'C': 108.5, 'Q': 143.8, 'E': 138.4, 'G': 60.1,\
               'H': 153.2, 'I': 166.7, 'L': 166.7, 'K': 168.6, 'M': 162.9, 'F': 189.9, 'P': 112.7, 'S': 89.0,\
               'T': 116.1, 'W': 227.8, 'Y': 193.6, 'V': 140.0}
hydrophobicities = {'A': 8.1, 'R': 10.5, 'N': 11.6, 'D': 13.0, 'C': 5.5, 'Q': 10.5, 'E': 12.3, 'G': 9.0,\
                    'H': 10.4, 'I': 5.2, 'L': 4.9, 'K': 11.3, 'M': 5.7, 'F': 5.2, 'P': 8.0, 'S': 9.2, 'T': 8.6,\
                    'W': 5.4, 'Y': 6.2, 'V': 5.9}
cmd.set('dot_solvent', 1)
cmd.set('dot_density', 3)
sasa = cmd.get_area(f'resi {resi} and chain {chain}')
rsa = sasa / maxASA[from_residue]
offset = -1.6
rsa_factor = +2.2
dV_factor = +0.59
dH_factor = -0.29
rsa_dV_factor = -0.54
rsa_dH_factor = +0.49
# calculate
dV = (volumes[to_residue] - volumes[from_residue]) / 100
dH = (hydrophobicities[to_residue] - hydrophobicities[from_residue])
simba = offset + rsa_factor * rsa + dV_factor * dV + dH_factor * dH + rsa_dV_factor * rsa * dV + rsa_dH_factor * rsa * dH
print -simba  # negative stabilised, positive destabilised

Human reasoning
However, the major problem with all the tools mentioned is that human reasoning is often required to read the literature and make a solid hypothesis.
A resource I cannot stress enough, for human mutations, its utility is gnomAD, which contains phenotypically neutral mutations: for example, one really ought to not claim a homozygous missense is pathogenic because it is causing a loss of function by destabilising the protein, yet there are homozygous truncations roaming around in the human population.
The paper describing the crystal structure used and those of the homologues may point out key features which would missed by automated servers, such as secondary-structure elements that are speculated to perform certain roles.
In fact, without doing any calculations seeing in PyMOL if a mutation causes a clash in the core of a protein, alters the active site, or changes a residue to a smaller one in a region known to alter conformation upon a given condition, is a lot more informative even if not quantitative...

You could run a molecular dynamics simulation for example with Gromacs or Charm. The problem with pymol is that the structure is static and you don't know how much tension that point mutation creates. First of all you would run an energy minimization to remove bad non-physical overlaps between atoms, then add solvent molecules, and then simulate your system first with position restraints for equilibration and then without for 1ns to 100ns+.
Even then you cannot be sure that what comes out is the 'final' result as sometimes an induced fold change can take much longer. It depends on the system and the energy barriers included and the time scales that you are interested in. If the protein is large it may have many different physiologically relevant conformations and the point mutation might only change the equilibrium distribution of those states or totally prevent a transition.
In case you have more than one crystal structure (in different conformations) you may introduce these point mutations in all of them and simulate them for a while. But it is not that trivial to setup such a simulation you should contact an expert. There are many things that can go wrong. Other software you could use are Modeller and Rosetta, if I am not mistaken.

You could give SAAPpred a try.
This tool uses machine learning (in the form of random forest) to predict the local effects of a single point mutation on the protein structure.
I did my PhD in the group and I used SAAPpred to find the most (predicted) damaging mutations of G6PD, to then ran extensive MD simulations to evaluate the final stability of those mutants.
Keep in mind that the tools doesn't return a structure, but just the possible/likely effects of that substitution.

Cool this looks good.

then ran extensive MD simulations to evaluate the final stability of those mutants..... Did You try to mutare the actual protein too ? Any interesting validation result ?


text:
GTF file format causing GeneMark-ES error
I am predicting genes using GeneMark-ES, I managed to generate a gtf file, but when I tried to extract the nt/protein sequence using get_sequence_from_GTF.pl, I am running into an error in the format of the gtf file, which I find puzzing as it was the gtf file was generated using gmes_petap.pl (an script in the GeneMark-ES package).
Here is the command I used and the error message which also prints the first line of the genemark.gtf file:
$ get_sequence_from_GTF.pl genemark.gtf DT_CAJQTU010000037_1.fasta

error in file format GTF: CAJQTU010000037.1 Desoria tigrina genome assembly, contig: utg000040l, whole genome shotgun sequence     GeneMark.hmm3   gene    2099    4039    .       +       .       gene_id "1_g"; transcript_id "1_t";

I think I have found the line in get_sequence_from_GTF.pl that is causing the error, but not being a perl monk, I am having troubles to read figure out what is wrong with my gtf file and more importantly figure out how to fix it.
...
if ( /^(\S+)\t\S+\t(\S+)\t(\d+)\t(\d+)\t.\t([\+-])\t([012.])\tgene_id \"(\S+)\"/ )
    {
    ... (some stuff)
    }
else
{
  if ( /^\s+$/ ) { next };
  if ( /^#/ ) { next };
  die( "error in file format GTF: $_\n" );
}
...

I'm afraid this is a clear bug in the tool you are using. The regular expression requires that all fields in the GTF file contain no whitespace. The \S+ means "one or more non-whitespace characters" but your "source" field, the second field in the GTF file, contains spaces: Desoria tigrina genome assembly, contig: utg000040l, whole genome shotgun sequence.
The thing is, that is a perfectly valid field for a GTF file! GTF fields can contain spaces, they just can't contain tabs. I tried to find a repository for the software to suggest a fix, but found only an empty one. In any case, try changing the regex you show to this:
#if ( /^(\S+)\t\S+\t(\S+)\t(\d+)\t(\d+)\t.\t([\+-])\t([012.])\tgene_id \"(\S+)\"/ )
if ( /^([^\t]+)\t[^\t]+\t([^\t]+)\t(\d+)\t(\d+)\t.\t([\+-])\t([012.])\tgene_id \"(\S+)\"/ )
    {
    ... (some stuff)
    }

I can't guarantee that this will work as expected, please test! The details will depend on what they do with the captured groups (the some stuff), but at least that regex will match and I think I have kept the numbering of the capture groups.

I agree with this 100%. I read the regex and didn't spot the reply. There are quite a few problems, the use of a bare "." just means 'any character' rather '.' which in this example would fit (albeit it might be intentional). The final section '; transcript_id "1_t";' isn't captured. However, the key problem is \S which will trip over the space characters in "CAJQTU010000037.1 Desoria tigrina genome assembly, contig: utg000040l, whole genome shotgun sequence". The use of '[^\t]', i.e. not a tab character will work. Personally I would split on the tab character.

Just to clarify this "\t.\t" means any single character between two tabs. However, this "[012.]" might mean a \. i.e. just a full-stop, I can't remember. If it doesn't [.] is exactly the same as [012.] or even [a-z|0-9] and beyond. Anyway this is the accepted answer.

You might use agat_sp_extract_sequences.pl from AGAT instead.
To extract and translate the coding regions:
agat_sp_extract_sequences.pl -g infile.gxf -f infile.fasta -p
see here for all possibilities with that script: https://agat.readthedocs.io/en/latest/tools/agat_sp_extract_sequences.html
Otherwise with the modification suggested by @terdon, your script should work fine.


text:
Why are my Nextflow processes not executing in parallel?
I have written a Nextflow script with three process:

The first process takes a pair of fastq files and aligns with reference genome. The process writes the resulting SAM file into sam channel.
Second process takes input from the sam channel and creates a BAM file from it, and writes it into bam channel.
Final process reads from bam channel and sorts the BAM files.

Following is the entire code:
#!~/bin nextflow

reads_ch = Channel.fromFilePairs('raw/*_{1,2}.fastq.gz', flat: true)
ref_genome = file('reference/genome/index/ref_genome')

process mapToRef {
    memory '20 GB'
    cpus 16
    input:
    tuple val(sample_id), file(read1), file(read2) from reads_ch

    output:
    tuple val(sample_id), file("${sample_id}.sam") into sam_ch

    script:
    """
    ~/biotools/bwa/bwa mem -t 60 ${ref_genome} ${read1} ${read2} > ${sample_id}.sam
    """
}

process samToBam {
    memory '20 GB'
    cpus 16
    input:
    tuple val(sample_id), file("${sample_id}.sam") from sam_ch

    output:
    tuple val(sample_id), file("${sample_id}.bam") into bam_ch

    script:
    """
    ~/biotools/SAMTOOLS/samtools-1.12/samtools view -S -b ${sample_id}.sam > ${sample_id}.bam
    """
}

process sortBam {
    memory '20 GB'
    cpus 16
    input:
    tuple val(sample_id), file("${sample_id}.bam") from bam_ch

    output:
    tuple val(sample_id), file("${sample_id}.sorted.bam") into sorted_bam_ch

    script:
    """
    ~/biotools/SAMTOOLS/samtools-1.12/samtools sort ${sample_id}.bam -o ${sample_id}.sorted.bam
    """
}

All three processes successfully execute their commands. However, they are executed sequentially. For example, the first process mapToRef runs sequentially for all read pairs, followed by samToBam, which too, runs sequentially after all mapToRefs have finished execution. I was under impression that parallelization would "work" out of the box, but I must be doing something wrong. I am a complete beginner in Nextflow and would greatly appreciate any help regarding this.
EDIT
Output of reads_ch.view():
[ERR2512393, /home/work/raw/ERR2512393_1.fastq.gz, /home/work/raw/ERR2512393_2.fastq.gz]
[ERR2512394, /home/work/raw/ERR2512394_1.fastq.gz, /home/work/raw/ERR2512394_2.fastq.gz]
[ERR2512391, /home/work/raw/ERR2512391_1.fastq.gz, /home/work/raw/ERR2512391_2.fastq.gz]
[ERR2512392, /home/work/raw/ERR2512392_1.fastq.gz, /home/work/raw/ERR2512392_2.fastq.gz]
[ERR2512390, /home/work/raw/ERR2512390_1.fastq.gz, /home/work/raw/ERR2512390_2.fastq.gz]

I have tested it on my end with the following script:
#!/usr/bin/env nextflow

reads_ch = Channel.fromFilePairs('/root/test_nextflow/*_{1,2}.fastq.gz', flat: true)

process mapToRef {
  memory '1 GB'
  echo true
  cpus 1

  input:
  tuple val(sample_id), file(r1), file(r2) from reads_ch

  output:
  tuple val(sample_id), file("${sample_id}.sam") into sam_ch

  script:
  """
  echo bwa mem ${r1} {r2}
  if [ ${sample_id} == "ERR2512393" ]; then sleep 10; fi
  touch ${sample_id}.sam
  """
}

process samToBam {
  memory '1 GB'
  cpus 1
  echo true

  input:
  tuple val(sample_id), file(sam) from sam_ch

  output:
  tuple val(sample_id), file("${sample_id}.bam") into bam_ch

  script:
  """
  echo samtools view -S -b ${sample_id}.sam
  touch ${sample_id}.bam
  """
}

process sortBam {
  memory '1 GB'
  cpus 1
  input:
  tuple val(sample_id), file(bam) from bam_ch

  output:
  tuple val(sample_id), file("${sample_id}.sorted.bam") into sorted_bam_ch

  script:
  """
  echo samtools sort ${sample_id}.bam
  touch ${sample_id}.sorted.bam
  """
}

It's producing the following result:

As you can see, the two processes are running in parallel (2 of 2)
It could be that yours are limited by the resources? If your PC doesn't have 40GB of memory and 32 CPUs, it will run them sequentially.
EDIT: By adding sleep 10 to mapToRef, we can clearly see both of them are running in parallel now:

EDIT 2: Added condition to "sleep" on one of the samples:

Final EDIT: be mindful of the resources you allocate to your processes. You could use ${task.cpus} and ${task.memory.toGiga()} to facilitate the configuration of your processes so that they are not different from what you have defined and what is actually used. See CPUs

Thanks @Law.
My execution screenshot also looks almost exactly the same.
Can you replace touch ... on your script with sleep 5 to emulate long running tasks and confirm if the processes are running parallelly? I have well over 20GB memory and 16 CPUs on this machine, so limited resource should not be an issue.

@PPrasai, check again my edit, i was doing it ;)

When the first process shows 1 of 2, does the second process start execution immediately? In my system, both reads from the first process need to be complete before the second process can start

Check my second edit. I have added a condition to mapToRef to "sleep" if one of the sample_id match one of the samples. As you can see, both of them runs in parallel first, since one of the samples doesn't match the condition, the other processes are executed afterward.

I seems to be working on your end...
Can you spot any issues with my code?

Try to lower the ressources specified. Also, in your Input channels, instead of specifying, for example, tuple val(sample_id), file("${sample_id}.sam") from sam_ch, you should rather do tuple val(sample_id), file(sam) from sam_ch and in the script, do ~/biotools/SAMTOOLS/samtools-1.12/samtools view -S -b ${sam} > ${sample_id}.bam

Also, I just noticed something... why is there 60 CPU defined in your mapToRef when you gave it 16 CPUs ?

Let us continue this discussion in chat.

Thanks @Law!
Reducing cpus to 4 on process and mapToRef, and memory to 4 GB make it all work parallelly! Thank you so much. Could you update your answer so that I can approve it

Well, I did say toward the end the following It could be that your are limited by the ressources? If you PC doesn't have 40GB of memory and 32 CPUs, he will run them sequentially. Also, you're welcome.


text:
If regulatory proteins decide the gene regulation, is it necessary for every cell to have a different look?
Transcription factors regulate the genes in different cells. That is, different cells have different sequence of genes that are transcribed (set on or off) by the RNA polymerase. This will help different types of cells to perform their own unique operation.
So, why do the cells have to look different physically, if the their function is determined by the on/off sequence?
NOTE: I am a neophyte in Bioinformatics. Any help on my building my conceptual knowledge is appreciated.

I’m voting to close this question because it belongs on biology SE.

There are many more contributors to the cell differentiation that the regulators. But one that you will encounter at an higher level is the chromatin status: if it is in the open status it allows transcripts to be on and being regulated (with all the further levels of regulations such as methylation, glycosylation etc). If the chromatin is off (the state is very condensed) very few transcript can be made from that region.
Check this wikipedia page here


text:
What is meant by transcriptional changes executed by the cell over a time period?
I read the following line in the research paper - The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells:

During differentiation, for example, each cell makes independent fate decisions by integrating a wide array of signals from other cells and executing a complex choreography of gene-regulatory changes. Thus individual cells can execute
the same sequence of transcriptional changes over highly varying
time scales.

Can anyone explain in simple terms what is meant by "same sequence of transcriptional changes"?
How can it be the same sequence if during differentiation, different types of cells are formed?

You are not mentioning, neither the paragraph, any specific time during development. The paragraph can be simply referring to stages in development when you already have groups of cells committing to the same transcriptional program (e.g. same part of an organ), as opposed to early stages where one single cell has to start a group of similar cells,

Not all cells will perform the same sequence of transcriptional changes. The paragraph you quote simply states that differentiation isn't a step change with all cells differentiating at the exact same moment. Rather, it's a more gradual change with cells each working on their own time scale.

That is right, embryogenesis is a good example in this case. Changes occur in the fetus when cells receive a mix of signal internal cellular factors as well as extracellular factors. However, the differentiation does not happen at the same time because the signals first. When some cells start to differentiate, they induce differentiation to the neighboring cells. Thus, a gradual change of differentiation.

Thanks @Devon Ryan for your time and expertise. Can you kindly explain what is meant by "sequence of transcriptional changes"?

A given order (sequence) of changes in gene expression.


text:
Does the alignment need to be redone when using SAW method to test long branch attraction?
I am looking into different methods to test the long-branch attraction problem. In this chapter they talk about the SAW method - excluding one and the other taxa from the dataset and see if they change their position within the tree. From the description of the method, I was not sure, if the bunch of sequences with excluded samples should be realigned, or if I should simply exclude the rows in the alignment.
I tried to figure out it by looking at the Siddall and Whiting 1995 paper, which was referred to as the original description of this method in the chapter mentioned above. The best I found is

Using their alignment (http://biology.unm.edu/esloker/), when Polypodium is excluded, the most parsimonious tree groups the myxozoans as sister to the Bilateria (Fig. 1a).

which makes me think that we should just exclude the taxa from the final alignment, but as there is still a space for interpretation, I thought I'll ask.

In my opinion, since the aim of the SAW test and other leave-one-out approaches like it is to minimise the variance associated with any factor other than the taxon sampling, you should leave the alignment as-is and simply exclude the row of the associated taxon.

Of course, things on long branches can often be annoyingly divergent and induce long internal INDELs and other weird alignment patterns. So I suspect it depends on exactly how much the homology statements in the alignment would change if you reinferred it dropping the taxon. Easiest way is just to try both (with realign, without realign) and compare.

That's already quite close to an answer, thanks @NatWH. (trying to subtly encourage you to write one)

I agree more with @NatWH's second comment. the idea of the test is to exclude all sources of information from the branch. Especially if the concern is long branches, the clustering that is implicit in the multiple alignment are likely to become nontrivial sources of information. One suspects also that choice of aligner (e.g. structure of initial heuristic steps) is also likely to interact with this.


text:
Subset GTF file for specific genes
I have a GTF file (in 'test.gtf' file):
chr1    HAVANA  gene    10069   18412   .   +   .   gene_id "ENSG00000223972.4"; transcript_id "ENSG00000223972.4"; gene_type "pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1"; level 2; havana_gene "OTTHUMG00000000961.2";
chr2    HAVANA  gene    11869   14412   .   +   .   gene_id "ENSG0000023172.1"; transcript_id "ENSG000023172.1"; gene_type "pseudogene"; gene_status "KNOWN"; gene_name "PCNA"; transcript_type "pseudogene"; transcript_status "KNOWN"; transcript_name "PCNA"; level 2; havana_gene "OTTHUMG00000001931.2";
chr3    HAVANA  gene    15769   19412   .   +   .   gene_id "ENSG00000433972.1"; transcript_id "ENSG000023172.1"; gene_type "pseudogene"; gene_status "KNOWN"; gene_name "USP21"; transcript_type "pseudogene"; transcript_status "KNOWN"; transcript_name "USP21"; level 2; havana_gene "OTTHUMG00000002361.2";

I have loaded it as a pyranges object and I want to subset it so that it only contains gene from test_list:
import pyranges as pr
import pandas as pd

gtf_file = "test.gtf"
test_list = ["PCNA","USP21", "USP1"]

gtf = pr.read_gtf(gtf_file)

sub_gtf = gtf

sub_gtf.df = gtf.df[gtf.df.gene_name.isin(test_list)]

But the last line throws this error:
File "<ipython-input-24-af0b4b5c7576>", line 1, in <module>
    irongtf.df = irongtf.df[irongtf.df.gene_name.isin(test_list)]

  File "/home/niek/.local/lib/python3.8/site-packages/pyranges/pyranges.py", line 314, in __setattr__
    _setattr(self, column_name, column)

  File "/home/niek/.local/lib/python3.8/site-packages/pyranges/methods/attr.py", line 48, in _setattr
    df.insert(pos, column_name, column)

  File "/home/niek/.local/lib/python3.8/site-packages/pandas/core/frame.py", line 4413, in insert
    value = self._sanitize_column(value)

  File "/home/niek/.local/lib/python3.8/site-packages/pandas/core/frame.py", line 4504, in _sanitize_column
    com.require_length_match(value, self.index)

  File "/home/niek/.local/lib/python3.8/site-packages/pandas/core/common.py", line 531, in require_length_match
    raise ValueError(

ValueError: Length of values (2) does not match length of index (1)

What would be the right way to subset the gtf file so that it only contains the rows where gene_name contains the gene symbols from test_list?
Or is there another tool that I could use to do this?

It turned out to be a lot easier if I used grep to do this:
grep -w -f genes.txt gencode.v19.annotation.gtf > sub_set.gtf

genes.txt contains gene symbols on new lines.

the only risk you are taking with this approach is that some of the gene names might be sub-strings of other gene names and as a consequence, those would get pulled as well.

I thought that the -w flag would prevent that?

Ah, you are right, I overlooked -w. Should be fine, unless the genes names contain non-word constituent characters (like - or .).

Maybe a better way would be: grep -wf <(awk '{print "gene_name \"$0\""}' genes.txt) gencode.v19.annotation,gtf > subset.gtf. This will ensure that the strings are compared only to the gene_name tag in each feature.

Great suggestion Ram RS; I have tried it but it gives me an empty file?

I'm not sure what the problem could be. Try just the awk command on the first 1000 lines of the big GTF, then try the grep with one of the output lines from the awk command. Step by step debugging like this is the way to solve these problems.

@justinian482 have you ever edited the genes.txt file in Windows? Could there be Windows line endings?

@terdon no I am using Ubuntu 20.04.3

Following up on zorbax's answer, you could read in and filter the GTF file in this way, among others:
#!/usr/bin/env python

import gtfparse as gp

gtf_file = "test.gtf"
test_list = ["PCNA", "USP21", "USP1"]

df = gp.read_gtf(gtf_file)
subset = df[df['gene_name'].str.contains('|'.join(test_list))]
print(subset)

You could further cull this subset by other attributes or columns, e.g.:
genes_only = subset[subset["feature"] == "gene"]

etc.
The gtfparse library is available via PyPi: https://pypi.org/project/gtfparse/ or Bioconda: https://bioconda.github.io/recipes/gtfparse/README.html

How would I save the df as a proper GTF after subsetting? The separators of the final columns are different compared to the first ones (; vs \t)?

You could add import pandas as pd and then try df.to_csv(out_filename, sep='\t') to write out a tab-delimited file from various data frame columns. You'll probably need to first reduce the attribute columns to a single column (https://stackoverflow.com/questions/27298178/concatenate-strings-from-several-rows-using-pandas-groupby), and then remove those additional columns. Or create a new data frame from the non-attribute and attribute columns and export that. Perhaps try it out and if you get stuck, ask another question showing what you've tried.

If you have your gtf file like a DataFrame you can use:
 df[df['gene_name'].str.contains('|'.join(test_list))]


text:
How do you set the coverage in PacBio's Sequel II?
I am reading the Whole Genome Sequencing for de novo Assembly Best Practices

Use the Sequel II or IIe System and SMRT® Cell 8M to sequence to desired coverage depth for complexity
of genome
10- to 15-fold coverage per haplotype recommended

First of all Wikipedia shows different definitions for the term haplotype. Which one applies here?
Second while the recommendations state to use 10- to 15-fold coverage per haplotype I find myself confused about how I would know the haplotype beforehand.
From what I understand of haplotype-estimation I would need the sequenced data first.
How would I go about this when wanting to do a de novo assembly of a species that has no published assembly of any close relatives? That means the only information I would have at that point is the C-value of the species to be sequenced.
Or better, how does the Sequel II figure out the number of haplotypes for achieving the desired fold-coverage?
I hope I'm in the right SE.

They just want to know the ploidy of the sample. In a diploid sample, e.g. humans, there are two haplotypes. This is mostly important for haplotype resolution (i.e. get the two full genomes in a diploid).
Note that elsewhere in the brochure that you link it writes:

Fast and efficient assembly of even
the largest genomes with haplotype
resolution of complex polyploids

If you are working with a diploid species, that translates to 20-30X coverage. That seems low to me for the older tech of CLR/Sequel sequencing, but they are the experts.
Different recommendations posted here by a service provider suggest 40-50X for structural variant analysis, which seems more reasonable.

There are still situations where I might not know the ploidity due to a lack of data. At least for now I was not able to find any data on either phylum: Cnidaria or  the class: Scyphozoa or the order: Rhizostomeae. Or in the case of Sequoia sempervirens how did they estimate the ploidity in advance?  Thanks for the insights!

@ilamengl I think that there are some methods you can look into: 1) karyotype analysis, 2) computational tools, here is one example: https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-018-2128-z

Additionally, there is a ploidy database which does appear to have cnidarians: https://cromanpa94.github.io/ACC/

Unless you want to reconstruct all haplotypes, you can sequence to ~30X regardless of ploidy.

I would agree for recent autopolyploids, but probably will be less true for polyploids with divergent haplotypes (e.g. allopolyploids).

@KamilSJaron There is not a good way to assemble those genomes for the moment and therefore few know what coverage is needed anyway. 30X would still be a good starting point.


text:
Why does my volcano plot show genes at the extremities of the graph?
I am confused about the genes at the top of the volcano plot. There are a lot of genes at the top, what's wrong with it? If that's normal what's the interpretation of it?
I used volcano plot tools from  galaxy with  p-value=0.05.

I suggest you give some background information about this plot as well as part of the data and the code that generates this plot. Asking this question without providing more information would be hard for us to help you. At first, I would say that you should set the Y axis to a higher value.

thank u @Law  for your comment, for the data I have three replicate of each condition, two conditions:  treated and control, in a total: 12 samples in fastq format. i used this pipeline :  Trimmomatic -> Kallisto ->  Deseq -> volcano plot  in galaxy

paired-end reads *

Either your graph or your data is capped at 300.  Probably the p-values themselves; whatever calculated those p-values literally can't represent a number tinier than 1x10^-300

thank u @swbarnes2  After trimming the quality of the data is good, there are no adapters or overexpressed sequences. for the alignment, I used Kallisto, then deseq for differential analysis.
I tried to augment the max of the axis of the plot,  but the top is 300. so us result we can't accepte this result ? The analysis will have to be redone ?

If the top is 300, you have to look for another galaxy website where the configuration of the volcano plot is more flexible. Or you can plot it manually with R. There are many websites that teaches you how to do a volcano plot. But what is intriguing to me is having a -log10 of more than 300... this seems excessive to me but it all depends on the research you are doing.


text:
blastp option to compare with a personal database
I'm trying to compare proteins generated with prokka (from a genome) with blastp using a virulence factor database VFDB, and I want to use some parameters obtained from this paper:
blastp option, top hit, ≥50% identity, ≥50% alignment length, E-value of < 10 - 5
How can I specify those parameters on the command line?
# create DB
makeblastdb -in VFDB.fas -dbtype prot

# run blast
blastp -query my.faa -db VFDB.fas -out results.txt .....

any suggestion?

You can use these options:
blastp -query my.faa -db VFDB.fas \
       -perc_identity 50 -outfmt 6 \
       -evalue 10e-5 -out results.txt

Then select the 5th column, qlen, and the 13th columns, slen, and get the percentage higher than 50.
awk  -F'\t' '{ if ($NF >= 50) printf "%s\t%.2f\n", $0 , ($2/$5)*100 }' result.txt

Check the documentation for the full -outfmt column names with blast -help

Thanks so much for the help. one more questions. I use the outfmt 6, and it have the next columns: qseqid, sseqid, pident, length, mismatch, gapopen, qstart, qend, sstart, send, evalue, bitscore. in this case which columns do I have to use to determinate the percent of alignment length? , thanks so much

you can configure outfmt 6 to produce a custom format, qlen and slen included, check the blastp -help


text:
What is the best QC to do on imputed UK Biobank data?
I am receiving imputed data from UK Biobank to conduct a GWAS on.  Previously I have carried out GWAS on genotype data, which I have QC'd for missingness per individual and per SNP, sex discrepancy, MAF filter < 0.05, heterozygosity rate, cryptic relatedness and population stratification.
However, from reading several papers about QC of imputed data, I do not carry out any of these QC checks on the imputed data.  What QC procedures should I expect to carry out on UKBB data? Just to filter by info score?
Or would it be better to carry out the initial QC checks on the genotype data and then impute myself?

In general my advice for stuff which is as subjective as QCing is to find a recent paper in a high profile journal (Nature Genetics perhaps) which has done a GWAS on imputed genotypes and see what they did in terms of QCing. There would be no need to impute yourself though.

A lot of effort went into the imputation of UK biobank by people who have been developing tools for imputation for many years. I don't think you'll be able to do better unless you are such a methods developer yourself or have access to a really nice secret reference panel. disclaimer: I developed tools that were likely used in the imputation of the uk biobank data that you are working with.

Thanks both for the comments, they were very reassuring!


text:
Does the Protein Data Bank contain the estimated distances obtained using NMR spectroscopy?
From the Wikipedia entry on the Protein Data Bank (PDB) (emphasis mine):

Most structures are determined by X-ray diffraction, but about 10% of structures are determined by protein NMR. When using X-ray diffraction, approximations of the coordinates of the atoms of the protein are obtained, whereas using NMR, the distance between pairs of atoms of the protein is estimated. The final conformation of the protein is obtained from NMR by solving a distance geometry problem.

I am currently working on an algorithm to solve distance geometry problems in the context of molecular conformations. I would prefer to test this algorithm on real data, so does the PDB contain information about the distances between pairs of atoms that were estimated from NMR spectroscopy measurements?
I checked the PDB file for the 2KB7 protein, but I only found the xyz coordinates of the atoms obtained after solving the distance geometry problem. Is it possible to get access to the estimated distances that were used to generate these conformations?

The data can be found here:
https://bmrb.io/search/instant.php?term=2KB7
There is no link directly in the PBD as it is a one-to-many relationship.
Clicking on the first value (a shifts dataset) you get somewhere in the middle a link to the actual dataset, which is a CIF like file with the following table:
                              vvvvvvv
    1 . 1 1  2  2 MET HA H  1   4.070 . . 1 . . . A  2 MET HA . 18256 1
    2 . 1 1  2  2 MET C  C 13 176.670 . . 1 . . . A  2 MET C  . 18256 1
                              ^^^^^^^

Each row has the spectral shift in ppm (chevrons added) of an assigned hydrogen or carbon atom. There are also anisotropic chemical shifts and dipolar couplings datasets etc. etc. Each with a particular piece of data...
Chemical structures are "elucidated" (=solved) in different ways depending on the type of experiment.
Say you wanted only the distanced calculated by Nuclear Overhauser effect (r_NOE), you would need to find a NOE NMR experiment.

Thanks a lot for your answer. I'm relatively new to the NMR literature. How are spectral/chemical shifts related to the distances between pairs of atoms in a protein?

Like the solutions to the phase problem in crystallography, the elucidation of the distances depends on the method used and differs between methods. However, the r_NOE mentioned above is a common technique. The elucidation of structures from a simple chemical shifts is not possible in a simple way without a myriad hacks.


text:
How to extract all sequences mapped to a transcript from Kallisto output
This question has also been asked on Biostars
I ran Kallisto with the --pseudobam option. How do I extract all the short reads that are mapped to a single transcript (e.g. ENST00000367969.8)?  As a person without any previous SAM/BAM experience, I tried the following things without success. Really appreciate some help here.
I first tried:
samtools view -S pseudoalignments.bam ENST00000367969.8

but got the following error:
[E::idx_find_and_load] Could not retrieve index file for 'pseudoalignments.bam'
[main_samview] random alignment retrieval only works for indexed BAM or CRAM files.

I tried to index the file using:
samtools index pseudoalignments.bam

but get the following error:
[E::hts_idx_push] NO_COOR reads not in a single block at the end 79729 -1
[E::sam_index] Read 'A00742:184:HLTF2DSXY:4:1101:4417:1000' with ref_name='ENST00000556431.1', ref_length=4142, flags=83, pos=2098 cannot be indexed
samtools index: failed to create index for "pseudoalignments.bam"

Sort the BAM file again with samtools sort and then run samtools index on the result. That should fix the issue.

What Devon said. Alternatively, you could also just convert the bam to sam and then look for lines where the 3rd field (the reference the sequence is aligned to) is what you want:
samtools view pseudoalignments.bam | awk '$3=="ENST00000367969.8"'

But really, it's far better to sort and index the file as Devon suggested so  you can use dedicated tools in future.


text:
Reproducible PyTorch Model
I'm using PyTorch (1.7.1), PyTorch Geometric (1.6.3), NVIDIA Cuda (11.2).
I need to make a neural network reproducible for a competition.
However, when I try:
device = torch.device('cuda:0')
rand = 123
torch.manual_seed(rand)
torch.cuda.manual_seed(rand)
torch.cuda.manual_seed_all(rand)
torch.backends.cudnn.deterministic = True
torch.backends.cudnn.benchmark = False
np.random.seed(rand)
random.seed(rand)

the results are different each time I run the code. How can I fix it?
(For a reference, it always comes out the same in device = torch.device('cpu').)

There is this nice function from Kaggle:
def seed_everything(seed=1234):
    random.seed(seed)
    os.environ['PYTHONHASHSEED'] = str(seed)
    np.random.seed(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed(seed)
    torch.backends.cudnn.deterministic = True

seed_everything()

There are several version of this function for torch, tensorflow and combinations. This specific one is from here.
Run this before you run your model. But not only the fitting function. It needs to be before the validation split function as well.


text:
Help to understand the code for dipeptide composition calculation (in python)
Dipeptide composition of a protein sequence is the number of times a particular dipeptide (e.g. Arginine-Histidine) occurs in a sequence divided by the total number of dipeptides in the sequence (which is the length of the sequence - 1)
I have found the following code for this:
import re

def DPC(fastas):
    AA = 'ACDEFGHIKLMNPQRSTVWY'
    encodings = []
    diPeptides = [aa1 + aa2 for aa1 in AA for aa2 in AA]
    header = ['#'] + diPeptides
    encodings.append(header)

    AADict = {}
    for i in range(len(AA)):
        AADict[AA[i]] = i

    for i in fastas:
        name, sequence = i[0], re.sub('-', '', i[1])
        code = [name]
        tmpCode = [0] * 400
        for j in range(len(sequence) - 2 + 1):
            tmpCode[AADict[sequence[j]] * 20 + AADict[sequence[j+1]]] = tmpCode[AADict[sequence[j]] * 20 + AADict[sequence[j+1]]] +1
        if sum(tmpCode) != 0:
            tmpCode = [i/sum(tmpCode) for i in tmpCode]
        code = code + tmpCode
        encodings.append(code)
    return encodings

There is another file (main) from where the arguments fastas come from which is basically a list of strings where every string is a protein sequence. I am finding it really difficult to understand the nested loops that were used:
    for i in fastas:
        name, sequence = i[0], re.sub('-', '', i[1])
        code = [name]
        tmpCode = [0] * 400
        for j in range(len(sequence) - 2 + 1):
            tmpCode[AADict[sequence[j]] * 20 + AADict[sequence[j+1]]] = tmpCode[AADict[sequence[j]] * 20 + AADict[sequence[j+1]]] +1
        if sum(tmpCode) != 0:
            tmpCode = [i/sum(tmpCode) for i in tmpCode]
        code = code + tmpCode
        encodings.append(code)

I have spent 4-5 hours on this but I am still finding it really difficult to follow. I would really appreciate if anyone can explain the steps involved in these loops. Thank you!
Edit: The github link to the files are :

DPC.py
iFeature.py

The github link to the python toolkit: iFeature

One thing that may help when trying to understand other peoples code is to add print statements and possibly break statements at various parts of the code.  You can print out exactly what is happening in the loop; for example, in the nested loop print j, print the tmpCode index, print the value at the tmp code etc.

@user438383, deleted

A few additional hints that may help clarify the algorithm:
The 20 magic number is because there are twenty amino acids.  The 400 magic number is because there are 20 * 20 potential dipeptides.
These peptides are then hashed by converting their letter to a number in the range [0, 20). You can then multiple one of them by 20 and add the other to get a number between [0, 400).  In the case of this function the first amino acid "code" is multiplied by 20. This number can be converted back to the dipeptide sequence when outputting the data.  The translation table they use is generated here diPeptides = [aa1 + aa2 for aa1 in AA for aa2 in AA].
The nested loops are to loop through each entry and then each sequence to get every dipeptide sequence, which can be done in linear time.
Just as a quick note on hashing DNA and peptide sequences.  You can generally convert them to an array index by using increasing powers of the alphabet size multiplied by the letter index of the alphabet.  For example, for DNA you can convert A=0, C=1, G=2, T=3, and you can hash fourmers by using 4^4 * a[i] + 4^3 * a[i + 1] + 4^2 * a[i + 2] + 4^1 * a[i + 3]. In code you may see this performed with bitshifts for DNA since the alphabet size is a power of 2.
The advantage of this hashing scheme is. thateach sequence is guaranteed to have a unique index.  The problem with this method is that as if you wish to use longer k-mers the memory usage will grow exponentially.  For example, for 16-mers you need 4 billion combinations, which is the size of a 32 bit unsigned integer.

Thank you so much. I am trying to follow what your suggestion. I am getting all sorts of error messages. I have very little programming skills. I only learnt to work with python on jupyter notebook. This toolkit (or any toolkit, I presume). works across multiple files at a time I presume. I am completely lost.

I have understood the concept of hashing i.e. AADict={'A': 0,  'C': 1,  'D': 2, 'E': 3, 'F': 4, 'G': 5, 'H': 6, 'I': 7, 'K': 8, 'L': 9, 'M': 10, 'N': 11, 'P': 12, 'Q': 13, 'R': 14, 'S': 15, 'T': 16, 'V': 17, 'W': 18, 'Y': 19}

I understood multiplying one of them by 20 and add the other to get a number between [0, 400). But still these 2 line if sum(tmpCode) != 0:             tmpCode = [i/sum(tmpCode) for i in tmpCode] are really bothering me

They are checking if the sum is larger than zero.  Then they are normalizing the totals to sum to 1.0; ie. calculating the frequency of each dipeptide.

It's working now. Thank you so much. I really appreciate


text:
Reason for linear effects dominating the difference in a quantitative trait between individual genomes
I'm an interested layman reading the paper https://arxiv.org/pdf/1408.3421.pdf (Hsu, Stephen, On the genetic architecture of intelligence and other quantitative traits, 2014). I'm trying to understand Hsu's claim (my possibly incomplete/faulty interpretation; for his words see pages 13 and 16):

Alleles affecting a quantitative phenotype trait that has been strongly selected for (such as height or intelligence for humans) will have their frequencies pushed close to 0 or 1 by evolution.
If there are two possible alleles at each of $N$ sites contributing to the quantitative trait, say with each minor allele having frequency $p$, then linear contributions of minor alleles would number $O(p)$ whereas non-linear effects would number $O(p^2)$. By the previous point, $p$ is small, so the linear terms dominate the total contribution.

Two things confuse me with this:

Might not the larger number of non-linear interactions (at least $\binom N2)$ vs $N$) compensate for the extra factor $p$?
Shouldn't we expect alleles only involved in non-linear contributions to experience less selective pressure, hence making the $p$ for such minor alleles larger than for those contributing linearly, again compensating for the differing orders in $p$?

Hoo boy. This is complicated. First, Hsu is a somewhat controversial character, many (most?) geneticists think of him as a scientific racist; his own students are trying to get him fired: https://statenews.com/article/2020/06/michigan-state-geu-calls-to-remove-vp-research-administrator-stephen-hsu?ct=content_open&cv=cbox_latest. That said, additive vs. epistatic effects is an active and complex debate. Your points 1 and 2 are reasonable and argued about frequently. Review articles arguing each way: https://pubmed.ncbi.nlm.nih.gov/24161664/, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4639398/

Thanks, those links are helpful!  I guess the argument that intelligence/height heredity being twice as high for monozygotic twins as for dizygtotic twins implies linearity seems convincing enough to me, although I don't get the more general theoretical argument put forward to explain why this should be the case.


text:
How to plot the gap distribution of contigs wrt to reference genome
I have a contigs file that I generated using Minia and then I have used Minimap2 to map these contigs again to the reference genome. Now I want to plot the gap distribution i.e. gap distance of the set of contigs with respect to the positions in the reference genome. For example: something like , if there are m contigs obtained from the Minia, then X axis will show 0 to m-1 gap points (in an ideal case where there is no overlap between the contigs and the contigs are aligned back to back without space).
The file from minimap2 looks like this:
0   17679   0   17679   +   NZ_BCML01000091.1   39426   18736   36415   17679   17679   60  NM:i:0  ms:i:35358  AS:i:35358  nn:i:0  tp:A:P  cm:i:3274   s1:i:17675  s2:i:49 de:f:0  rl:i:75 cg:Z:17679M
1   16011   0   16011   +   NZ_BCML01000021.1   105137  86      16097   16011   16011   60  NM:i:0  ms:i:32022  AS:i:32022  nn:i:0  tp:A:P  cm:i:3008   s1:i:15996  s2:i:0  de:f:0  rl:i:15 cg:Z:16011M
2   2779    0   2779    -   NZ_BCML01000048.1   72397   14902   17681   2779    2779    60  NM:i:0  ms:i:5558   AS:i:5558   nn:i:0  tp:A:P  cm:i:519    s1:i:2751   s2:i:0  de:f:0  rl:i:45 cg:Z:2779M
4   177     0   177     +   NZ_BCML01000150.1   16425   15710   15887   177     177     60  NM:i:0  ms:i:354    AS:i:354    nn:i:0  tp:A:P  cm:i:30     s1:i:160    s2:i:0  de:f:0  rl:i:15 cg:Z:177M

where the 1st column is the contig ID and 8th and 9th column show the start and end position of the contig on the reference genome.
I have used the samtools for obtaining a length distribution plot, but not understanding how to do a gap distribution plot. I can use awk on linux or python or R. Anything that is going to work. Thank you for all the help.

a better solution
I wasn't too happy with the previous answer that I gave, so I was fiddling around, and I realized that bedtools genomecov will give you this for free with the following workflow:
minimap2 -ax asm10 ref.fasta contigs.fasta | \
samtools sort --output_fmt BAM > contigs_to_ref.bam
# -bga means output in bedgraph format, including all positions
bedtools genomecov -ibam contigs_to_ref.bam -bga > contig_cov.bg

If you look at the resulting file, it will report all regions of identical coverage:
head contig_cov.bg
# prints the following:
# reference chr/contig, start position, end position, coverage
NODE_1_length_442924_cov_598.199    0   8151    0
NODE_1_length_442924_cov_598.199    8151    8392    1
NODE_1_length_442924_cov_598.199    8392    20203   0
NODE_1_length_442924_cov_598.199    20203   20463   1
NODE_1_length_442924_cov_598.199    20463   44953   0
NODE_1_length_442924_cov_598.199    44953   45216   1
NODE_1_length_442924_cov_598.199    45216   52633   0
NODE_1_length_442924_cov_598.199    52633   52899   1
NODE_1_length_442924_cov_598.199    52899   57698   0
NODE_1_length_442924_cov_598.199    57698   58207   1

You can simply filter this file for all regions wherein the fourth column is equal to zero, and take the first 3 columns, and that is a list of gaps in the assembly relative to the reference:
awk '$4 == "0" { print $0 }' contig_cov.bg | cut -f1-3
# prints
NODE_1_length_442924_cov_598.199    0   8151
NODE_1_length_442924_cov_598.199    8392    20203
NODE_1_length_442924_cov_598.199    20463   44953
NODE_1_length_442924_cov_598.199    45216   52633
NODE_1_length_442924_cov_598.199    52899   57698
NODE_1_length_442924_cov_598.199    58207   73398
NODE_1_length_442924_cov_598.199    73733   78601
NODE_1_length_442924_cov_598.199    78959   80147
NODE_1_length_442924_cov_598.199    80419   86967
NODE_1_length_442924_cov_598.199    87223   93872

You could even do this as a one-liner probably though that's getting a little convoluted.
Original answer
It's a bit circuitous, but one way to do this would be to follow this answer's suggestion and generate a pileup for the alignment of your contigs to the reference.
You can get minimap2 to give SAM/BAM output with the -a flag. You then use the SAM/BAM to generate a pileup or a coverage file, using samtools mpileup or samtools depth -a. You can then use a heuristic to find contiguous regions with zero coverage of contigs. You could for example loop through a depth file, here is some python code that gives the basic idea:
#!/usr/bin/env python

# minimap2 -ax asm10 ref.fasta contigs.fasta | \
# samtools sort --output_fmt BAM > contigs_to_ref.bam
# samtools depth contigs_to_ref.bam > depth.txt
depth_file = "depth.txt"

with open(depth_file) as dep_handle:
  gap_sizes = list()
  start = None
  pos = None
  for line in dep_handle:
    # needs also logic to handle multiple chromosomes in reference
    # but I'm lazy
    fields = line.split()
    if int(fields[2]) == 0:
      pos = int(fields[2])
      if start is None:
        start = int(fields[1])

    elif fields[2] >0:
      if start is None:
        continue
      else:
        gap_sizes.append(pos - start)
        start = None

print(gap_sizes)

Probably there is a simpler way to do this but I don't know it.


text:
Fastest way to count where two bam files are homozygote reference (inverted variant calling)
Maybe I'm having a brain-freeze so excuse me if this is a waste of everyone's time...
I am working with an organism with ~700Mb genome.
I have bam files of two individuals that are whole genome sequenced with the average depth of ~150.
I would like to count the occurrences when both of the individuals are homozygote for the reference allele. Preferably without any hassle or space, I have no need to create a vcf file for the task or save the information in any form. Something like 'gatk or mpileup... argument that only retains sites when all are homozygote reference... delete tmp file...'. So the question is: what would be the "cheapest"/fastest way to count where these individuals are homozygote reference?
Come to think of it... I'm asking for a line count value on an "inverted variant file"/"gVCF minus variation"

There is no genome reference for your organism ? What you want to do is very simple, you want to do variant calling using the BAM files (there are a number of ways for that) and then count the sites with variation, then you just do (number_of_all_sites_in_genome minus number_of_variant_sites)

I was thinking of all "callable" sites according to my bams (that meet some filtering criteria). Not just reference size minus variant sites. Maybe I will just do some genotype file and parse through it. I just thought I was forgetting some obvious command in the known tools that would be fast and easy.

Like bcftools call --skip-variants

gatk has an option to output all sites and not only variants, you could do your calculation on that VCF and you will be able to filter reference sites based on reference allele quality, on the other hand you will have to deal with a massive VCF and slow computation,  as opposed to the size minus #_variants approach,

For now it's gvcf via bcftools and related filtering and then something like awk if substr$x,1 == substr$x,3 and $x==$y.... Heavy duty for minor reward, the antonym of judo.