text:
Remove Redundant Sequences from FASTA file with Biopython reducing memory footprint
I need help with my Biopython code, I am kind of new to Python and Biopython so please don't be harsh on me.
Everything started from here: Remove duplicated sequences in FASTA with Python
I answered copying code from here: Remove Redundant Sequences from FASTA file in Python but I am still bugged about the OP's first question: Remove duplicated sequences in FASTA with Python
So I decided to try to minimize the code footprint on RAM. I ended up copying from Biopython documentation Split large file and producing:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Nov  11 17:40:03 2020

@author: Pietro

"""

from memory_profiler import memory_usage

import memory_profiler

def batch_iterator(iterator, batch_size):
    """ from https://biopython.org/wiki/Split_large_file
    """
    entry = True  # Make sure we loop once
    while entry:
        batch = []
        while len(batch) < batch_size:
            try:
                entry = iterator.__next__()
            except StopIteration:
                entry = None
            if entry is None:
                # End of file
                break
            batch.append(entry)
        if batch:
            yield batch

# @profile
def f():

    from Bio import SeqIO

    import time

    import sys

    start = time.time()

    data_input_iter = SeqIO.parse("Newtextfile-03", "fasta")
    print('data_input_iter : ',sys.getsizeof(data_input_iter),'\n')

    c=0
    for i , batch in enumerate(batch_iterator(data_input_iter , 1)):
        record_iter = batch

        val = True
        c += 1
        if c == 1:
            with open("removed-double_Newtextfile-03", 'w+') as output_file:
                    SeqIO.write(record_iter, output_file, "fasta")
                    # print('copy   : ',(record_iter))
                    output_file.flush()
            output_file.close()

        else:

            SeqIO.write(record_iter, "temp_file", "fasta")

            record_iter = SeqIO.read("temp_file", "fasta")

            data_output_iter = SeqIO.parse("removed-double_Newtextfile-03", "fasta")

            d = 0

            try:
                for i , batch in enumerate(batch_iterator(data_output_iter  , 1)):
                        SeqIO.write(batch, "temp_file2", "fasta")

                        record_output = SeqIO.read("temp_file2", "fasta")

                        d += 1
                        if record_iter.seq != record_output.seq  :
                            # print('pass   : ',record_iter.name)
                            pass
                        else:
                            val = False
                            # print('false')
                            break

            except:
                print('ERROR !!!!!!!!!!!')
                quit()

            if val:
                        # print('copy   : ',record_iter.name)
                        with open("removed-double_Newtextfile-03", 'a+') as output_file:
                            output_file.write(str(record_iter.format("fasta")))
                        output_file.close()
            del data_output_iter

    end = time.time()

    print(f"Run time is {(end- start)}")

mem_usage = memory_usage(f)
# print('Memory usage (in chunks of .1 seconds): %s' % mem_usage)
print('Maximum memory usage: %s' % max(mem_usage))

# f()

but I have to major problems :
FIRST: I wasn't able to find a way to pass a SeqRecord object from a list to a file without saving it to a temp file first. That is without:
for i , batch in enumerate(batch_iterator(data_input_iter , 1)):
        record_iter = batch
......

SeqIO.write(record_iter, "temp_file", "fasta")

record_iter = SeqIO.read("temp_file", "fasta")

the line:
 if record_iter.seq != record_output.seq  :

will throw:
AttributeError: 'list' object has no attribute 'format'

AttributeError: 'list' object has no attribute 'seq'

you get first error and the if remove format from:
output_file.write(str(record_iter.format("fasta")))

and #quit() you get the second error
How should I cope with this ? I mean is there a way to bypass the save and read lines ?
SECOND : and more important. I am sure my code doesn't use a set and a list .... and that's why is extremely slow compared to the copied Remove Redundant Sequences from FASTA file in Python posted below too. But how does it use even more RAM that the fastest code ? I tried to minimize the memory footprint but wasn't able !! Am I using a bad way to try to calculate the amount of RAM used ? And does Biopythom SeqIO.parse object already is the optimal choice for , lets say "2/4 Gb size" fasta files on a laptop with just 4Gb of RAM (remember the original OP post  Remove duplicated sequences in FASTA with Python). Maybe my test file is not big enough ? Does anybody have a 2GB fasta file with a dozen of duplicates in it ?
Remove duplicated sequences in FASTA with Python memory usage:
size of parser :  48
memory seen :  32984
memory records :  2888
Run time is 0.0008047382036844889
Maximum memory usage: 55.359375

my code memory usage and speed:
data_input_iter :  48

Run time is 37.781240463256836
Maximum memory usage: 56.40234375

My fasta test set is: https://easyupload.io/z2yss1
Biopython method to remove duplicated from  Remove duplicated sequences in FASTA with Python:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
copied from https://bioinformatics.stackexchange.com/questions/13708/remove-redundant-sequences-from-fasta-file-in-python

uses set() instead of list for seen

should work faster for big files

"""

from memory_profiler import memory_usage

def f():

    import sys

    from Bio import SeqIO
    import time

    start = time.time()

    seen = set()
    records = []

    parser = SeqIO.parse("Newtextfile-03", "fasta")
    for record in parser:
        # print('running')
        if str(record.seq) not in seen:
            seen.add(str(record.seq))
            records.append(record)
            # print('copying')

    print('size of parser : ',sys.getsizeof(parser) )
    print('memory seen : ',sys.getsizeof(seen))
    print('memory records : ',sys.getsizeof(records))

    #writing to a fasta file
    SeqIO.write(records, "removed-double_Newtextfile-03-from-002", "fasta")
    end = time.time()

    print(f"Run time is {(end- start)/60}")

mem_usage = memory_usage(f)
# print('Memory usage (in chunks of .1 seconds): %s' % mem_usage)
print('Maximum memory usage: %s' % max(mem_usage))

# f()

I wrote the answer here and while I'm not clear on all aspects of this question, I can offer two suggestions to reduce the memory footprint considerably while maintaining high speed.
First, rather than appending the records to a list you could use a generator that generates elements on-the-fly rather than storing them in memory, e.g.
seen = set()

def yield_seqs():
    for record in SeqIO.parse("b4r2.fasta", "fasta"):
        if record.seq not in seen:
            seen.add(record.seq)
            yield record

SeqIO.write(yield_seqs(), "output.fasta", "fasta")

Second you could save even more memory, but with the risk of false positives, by using a bloom filter, e.g. using python-bloom-filter, just replace the set with seen = BloomFilter(max_elements=10000, error_rate=0.001)

thanks for the generator hint.

FIRST: I wasn't able to find a way to pass a SeqRecord object from a list to a file without saving it to a temp file first. That is without:
for i , batch in enumerate(batch_iterator(data_input_iter , 1)):
        record_iter = batch
......

SeqIO.write(record_iter, "temp_file", "fasta")

record_iter = SeqIO.read("temp_file", "fasta")

the line:
 if record_iter.seq != record_output.seq  :

will throw:

AttributeError: 'list' object has no attribute 'seq'

After spending some timde on it I realized that:
print('\n@@@@@@@@@@@@@@\n',type(record_iter),'\n',record_iter)

is returning:
<class 'list'>

but jus using :
for i , batch in enumerate(batch_iterator(data_input_iter , 1)):
        record_iter = batch[0]

now :
print('\n@@@@@@@@@@@@@@\n',type(record_iter),'\n',record_iter)

is returning:
<class 'Bio.SeqRecord.SeqRecord'>

no need to save and read the Bio.SeqRecord.SeqRecord' anymore
so this is just half the question answered so no need to save and read the record on file
but even if speed increases the memory usage goes up :
data_input_iter :  48

Run time is 1.442115306854248
Maximum memory usage: 110.09765625
```


text:
Bulk RNA-Seq Read Length Normalization across different samples
I have 20 samples out of which 14 are 100 bp in length and 6 are 150 bp. Is there a way to normalize the read length for cross-sample differential expression comparison? What would be the best way to approach this?
Note: The samples are NOT from reference/model species.

Are these paired-end reads (i.e. 2x100bp; 2x150bp), or single-end reads? That will make quite a big difference in terms of the mappability of the reads to transcripts.

The simplest thing to do is to trim the 150 bp fasts so that they are 100 long.  I don't think there is an easy way to correct for the fact that the 150 bp long reads will have a higher unambiguous rate of alignment and gene assignment than the 100 bp long reads.
If you have a mix of all the experimental conditions across all the lengths, you can include length as an element in the design, so that DESeq2 may be able to model the difference between the two batches.

I've done this by doing a RPkM-like transformation which I call VSTPk. Transformations that adjust for transcript length will correct a little bit for the noise associated with mapping short reads to transcripts of varying lengths:
https://bioinformatics.stackexchange.com/a/15218/73
Another, probably better, alternative is to use Salmon to map reads to transcripts, which applies its own mapping-based correction to read counts. Salmon produces a statistic that relates to the mappable length of a transcript, and normalises counts based on how many molecules the model predicts came from each transcript isoform.
Unfortunately all these normalisation techniques require good gene models, which are difficult to obtain for non-model organisms.

I think this is not what OP is asking. OP has different read lengths, that has nothing to do with gene- or transcript length which RPKM-ish transformations are acting on. The concern here is mappability bias, not the fact that longer transcript inherently produce larger counts. I would therefore go with the suggestion above to simply trim back reads to the same length, that is both simple and effective to eliminate this potential bias without actually applying any data transformation or normalization, plus is independent of gene models and reference annotations. Agree with use of salmon.

source: Stack Exchange

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

sort: votes


text:
Mrbayes stop at certain generation without error
I am using Mrbayes v 3.2.7 MPI version to infer phylogeny, and I set the ngeneration to 10,000,000, but it often stop at certain generation (usually far smaller than given generations). Even if I restart the analysis by invoking the append=yes , but it stops again.
Anyone have this problem? I would appreciate you for sharing your useful experience.
Best regards,
Yanpeng
Version info

Version:   3.2.7
Features:  SSE AVX FMA Beagle readline
Host type: x86_64-conda-linux-gnu (CPU: x86_64)
Compiler:  gnu 7.5.0

Command line:
nohup mpirun -np 8 mb-mpi concatenate_phylipID.nex < /dev/null > /dev/null &

Mrbayes block in nexus file:
begin mrbayes;

    outgroup t1;

    charset Subset1 = 1-373 1953-2032;
    charset Subset2 = 374-531;
    charset Subset3 = 532-677 1583-1618;
    charset Subset4 = 678-1442;
    charset Subset5 = 1443-1582 1619-1725;
    charset Subset6 = 1726-1802 1803-1952;
    charset Subset7 = 2033-2661;

    partition PartitionFinder = 7:Subset1, Subset2, Subset3, Subset4, Subset5, Subset6, Subset7;
    set partition=PartitionFinder;

    lset applyto=(1) nst=6 rates=invgamma;
    lset applyto=(2) nst=6 rates=invgamma;
    prset applyto=(2) statefreqpr=fixed(equal);
    lset applyto=(3) nst=6 rates=invgamma;
    lset applyto=(4) nst=6 rates=invgamma;
    lset applyto=(5) nst=6 rates=invgamma;
    lset applyto=(6) nst=6 rates=gamma;
    lset applyto=(7) nst=6 rates=invgamma;

    prset applyto=(all) ratepr=variable;
    unlink statefreq=(all) revmat=(all) shape=(all) pinvar=(all) tratio=(all);
    mcmc nchains = 4 ngen = 10000000 samplefreq = 10000 printfreq = 10000 diagnfreq = 10000;
    sump;
    sumt;
    quit;
end;

How are you deciding it has stopped prior to 10,000,000 generations? Is it by counting the number of samples?

Check the log file,  It tells you how many generations have been done.

Have you tried running it in a screen or similar to see if it is something weird with nohup?

I prefer running prpgroms background if the analysis will take a long time. I think you are right. Ever I run it in PC five days, it works. Thanks.

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


text:
How to grep with a list of genes from a gff the GO terms
I have a text file with a list of name of genes (alias names) like this (** added for emphasis):
Solyc00g019980
**Solyc01g005020**
..
.

And i have a gff file like this:
SL4.0ch01   maker_ITAG  gene    46201   59605   .   -   .   ID=gene:Solyc01g005020.3;Alias=**Solyc01g005020**;Name=Solyc01g005020.3;length=13404
SL4.0ch01   maker_ITAG  mRNA    46201   59605   .   -   .   ID=mRNA:Solyc01g005020.3.1;Parent=gene:Solyc01g005020.3;Name=Solyc01g005020.3.1;Note=Salt Overly Sensitive 1;_AED=0.25;_QI=241|1|1|1|0|0|23|311|1151;_eAED=0.25;**Ontology_term=GO:0006812,GO:0015299,GO:0016021,GO:0055085**
SL4.0ch01   maker_ITAG  three_prime_UTR 46201   46511   .   -   .   ID=three_prime_UTR:Solyc01g005020.3.1.0;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  exon    46201   46844   .   -   .   ID=exon:Solyc01g005020.3.1.1;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  CDS 46512   46844   .   -   0   ID=CDS:Solyc01g005020.3.1.1;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  exon    47281   47585   .   -   .   ID=exon:Solyc01g005020.3.1.2;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  CDS 47281   47585   .   -   2   ID=CDS:Solyc01g005020.3.1.2;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  exon    48276   48543   .   -   .   ID=exon:Solyc01g005020.3.1.3;Parent=mRNA:Solyc01g005020.3.1

I would like to grab only the ones i have in my list gene text file (Alias names) and the corresponding "Ontology_term" next to it.
Terminal is more convenient if there is a solution like grep or something.
Thanks in advance
edit:
i would like to get this output, as an example
Solyc01g005070  Ontology_term=GO:0005543
Solyc01g005230  Ontology_term=GO:0008168

i dont mind if this is included "Ontology_term=", i would prefer not to be included though
thanks again
grep -f list.txt data.gff doesn't work for me, because it grabs all the other as well.

Are you trying to grep the Ontology term from the next line in your example? If so, that makes things a little more complicated.

@JimHat your output example does not correspond to your input examples. Solyc01g005070 and Solyc01g005230 are not referred anywhere in your gff example. Thus, it's hard to relate exactly what you want exactly

@Jim Hat
I took a stab at this given the example files you included.
If the 'list.txt' file looks like this:
Solyc00g019980
Solyc01g005020

And the data.gff file looks like this:
SL4.0ch01   maker_ITAG  gene    46201   59605   .   -   .   ID=gene:Solyc01g005020.3;Alias=Solyc01g005020;Name=Solyc01g005020.3;length=13404
SL4.0ch01   maker_ITAG  mRNA    46201   59605   .   -   .   ID=mRNA:Solyc01g005020.3.1;Parent=gene:Solyc01g005020.3;Name=Solyc01g005020.3.1;Note=Salt Overly Sensitive 1;_AED=0.25;_QI=241|1|1|1|0|0|23|311|1151;_eAED=0.25;Ontology_term=GO:0006812,GO:0015299,GO:0016021,GO:0055085
SL4.0ch01   maker_ITAG  three_prime_UTR 46201   46511   .   -   .   ID=three_prime_UTR:Solyc01g005020.3.1.0;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  exon    46201   46844   .   -   .   ID=exon:Solyc01g005020.3.1.1;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  CDS 46512   46844   .   -   0   ID=CDS:Solyc01g005020.3.1.1;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  exon    47281   47585   .   -   .   ID=exon:Solyc01g005020.3.1.2;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  CDS 47281   47585   .   -   2   ID=CDS:Solyc01g005020.3.1.2;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  exon    48276   48543   .   -   .   ID=exon:Solyc01g005020.3.1.3;Parent=mRNA:Solyc01g005020.3.1

(Note I removed the **s, which I assumed were added for emphasis)
You can do the following:
aliasNames=$(grep -f list.txt data.gff | cut -d':' -f2 | cut -d';' -f1)
ontologyTerms=$(grep -f list.txt data.gff | cut -d';' -f 8 | cut -d'=' -f2)
paste <(echo "$aliasNames") <(echo "$ontologyTerms")

The first line filters the file based on your list and then extracts the alias names from the filtered GFF file. The second line filters the GFF and extracts the ontology terms (after the equal sign). And the third line joins the two as columns.
The result that is printed to the screen is:
Solyc01g005020.3
Solyc01g005020.3.1      GO:0006812,GO:0015299,GO:0016021,GO:0055085
Solyc01g005020.3.1.0
Solyc01g005020.3.1.1
Solyc01g005020.3.1.1
Solyc01g005020.3.1.2
Solyc01g005020.3.1.2
Solyc01g005020.3.1.3

Given the examples you included all of the lines of the GFF are included, but not all have gene ontology information.


text:
How do I keep a glob from copying a file in Cromwell / WDL
Using Crowell 52.  When I do a glob to get matching files, Cromwell places a copy of each matching file in the glob-* directory.
As these are multi-gigabyte Bam files, I'd really prefer for Cromwell to store soft links to the files, rather than full copies.
The code:
   output {
      Array[File] Outputs = glob("*.fastq.gz")
   }

Is it possible to make this happen?

It would help your case to show the code you're using.

@RamRS added above

As stated in https://github.com/broadinstitute/cromwell/issues/4395
The problem may be because you can't hardlink folders on your system for various reasons (i.e. filesystem doesn't allow hardlinking folders, only files,  or you are trying to link between different filesystems).
You can try overriding the default way glob works.
i.e. add this to your cromwell.config
glob-link-command = "ls -L GLOB_PATTERN 2> /dev/null | xargs -I ? ln -s ? GLOB_DIRECTORY"

Like this
backend
{
  default = "Slurm"
  providers
  {
    Htcondor {
      actor-factory = "cromwell.backend.impl.sfs.config.ConfigBackendLifecycleActorFactory"

      config
      {
        glob-link-command = "find GLOB_PATTERN 2> /dev/null | xargs -I ? ln -s ? GLOB_DIRECTORY"


text:
What is a PacBio "movie file"?
I came across references to "movie files" from PacBio sequencing in this paper:
https://www.jimmunol.org/content/204/12/3434
Specifically:

Movie files used to generate results presented in this article have been submitted to the National Center for Biotechnology Information BioProject (https://www.ncbi.nlm.nih.gov/bioproject/) under accession number PRJNA389440.

Looking at the SRA entries, it looks like it's just the raw reads as I expected.  But what is meant by "movie files"? (This is probably just my ignorance about PacBio as opposed to Illumina or others but I've never seen that term before and can't find much online.)

[disclaimer: not a PacBio person, I'm just vaguely aware of the technology]
Movie files are the PacBio raw data format, representing observed signal intensities in each sequencing well over the course of the sequencing run. The movie files can be used to re-call sequences with updated statistical models, which can improve accuracy, or test out hypotheses about how epigenetics alters the rate of nucleotide incorporation (i.e. methylation).
Illumina has a similar static version of this which is an image (or multiple images, depending on chemistry) of the flow cell after each synthesis cycle. Illumina's raw data format consumes tens of terabytes per run, and is of minimal value after base calling, so is almost always discarded after base calling.
Anything uploaded to NCBI SRA must contain FASTQ information. While the original uploaded files can be kept (which may include time/signal information), the files are converted into a FASTQ-like format for NCBI's internal representation. As such, there's a chance that movie files could be uploaded to SRA, but those files would also need to include the called FASTQ information (e.g. in a HDF5 container). In contrast to this, ENA does allow uploading raw data files without called information.

Got it!  In that paper's context do the authors probably mean they deposited reads in the SRA post-base-calling, then, rather than pre?  (I deposit FASTQ from my Illumina runs, not the raw images, for example.)  On the other hand I see NCBI's docs mention HDF5, so could it include signal intensity?


text:
How to extract column from a matrix matching the another file (sample file)
I have two files master file.txt and sample.txt
master.txt
Name    GTEX.1117F.3226.SM.5N9CT    GTEX.111FC.3126.SM.5GZZ2    GTEX.111FC.3326.SM.5GZYV    GTEX.1128S.2726.SM.5H12C    GTEX.1128S.2826.SM.5N9DI    GTEX.117XS.3026.SM.5N9CA
ENSG00000223972 1   0   0   1   0   0
ENSG00000227232 298 168 197 106 221 184
ENSG00000278267 0   1   1   0   0   0
ENSG00000243485 0   2   0   1   0   1
ENSG00000237613 0   0   0   1   0   1
ENSG00000268020 1   0   0   1   0   1
ENSG00000240361 1   3   0   4   2   1
ENSG00000186092 3   1   1   1   1   2

Sample file.txt
GTEX.1117F
GTEX.111FC
GTEX.111XS

Desired output
Name    GTEX.1117F.3226.SM.5N9CT    GTEX.111FC.3126.SM.5GZZ2    GTEX.111FC.3326.SM.5GZYV    GTEX.117XS.3026.SM.5N9CA
ENSG00000223972 1   0   0   0
ENSG00000227232 298 168 197 184
ENSG00000278267 0   1   1   0
ENSG00000243485 0   2   0   1
ENSG00000237613 0   0   0   1
ENSG00000268020 1   0   0   1
ENSG00000240361 1   3   0   1
ENSG00000186092 3   1   1   2
ENSG00000238009 2   2   1   1
ENSG00000233750 0   3   3   0
ENSG00000268903 103 44  76  24

I tried using the R-code form this link
Print specific columns in a matrix on the basis of sample id's in the header
However, I am not able to figure out how to extract columns matching only half part of ID (like for GTEX.111FC I want both GTEX.111FC.3126.SM.5GZZ2 and   GTEX.111FC.3326.SM.5GZYV)
Thank you

Assuming your "master" files looks like this:
df <- data.frame(Name = LETTERS[1:5],
                 Col1 = sample(1:10,5),
                 Col2 = sample(1:10,5),
                 Col11 = sample(1:10,5),
                 Col3 = sample(1:10,5))

  Name Col1 Col2 Col11 Col3
1    A    5    4     3   10
2    B    1    6     5    6
3    C    3    1     9    4
4    D    9    9     8    9
5    E    8   10     1    1

For a sequence of pattern you want to extract, you can use base r function such as grep:
match_pat = c("Col1","Col2")

m <- unlist(sapply(match_pat,function(x) grep(x,colnames(df))))

Col11 Col12  Col2
    2     4     3

And then, extract your columns of interest by doing:
> df[,m]
  Col1 Col11 Col2
1    5     3    4
2    1     5    6
3    3     9    1
4    9     8    9
5    8     1   10

In your example, it will be something like:
match_pat <- c("GTEX.1117F","GTEX.111FC","GTEX.111XS")
m <- unlist(sapply(match_pat,function(x) grep(x,colnames(df))))

master[,m]
```

Thank you for your answer. However, my actual sample list has more than 100 id's and matrix has 2000 rows. I want match_pat to read a list of id's? Is it possible?

@Megha where does your list of ID's come from ? a file ?

@Law: I have a GTEX brain sample file (which have id's in the header like  GTEX.111FC.3126.SM.5GZZ2. I want to extract the age-related samples(20-60 yrs) from this file. The age-related data file is a list of files which says all id's starting with GTEX.111FC are 20-30 year old and so on. The issue is the id's in the list is the half part of ID in the main master file making it difficult for me to extract it

And yes age id file is a separated file with the list of ID's

Hi, I did it with bash. I got the full id's from the partial match.
xargs -i grep ^{} matched.Brain.txt.2nd.txt < 2ndcol.20_30 > new.file

Glad you succeed to solve your issue with Bash. To answer your question, you can define match_pat to be a vector of the column containing id files like `match_pat <- as.vector(list_id_file$NameOfTheColumWithID). Then, you can extract column form your GTEX brain sample file.


text:
Denovo, Stacks: Getting an “ambiguous redirect” error
I tried running the following command
/public/home/WHC_BBT/Softwares/bin/bin/denovo_map.pl --samples $RawData -O $Popmap -o $fullM2n3 -M 2 -n 3 -m 3 -T 15 &> $logfile

This is the error I got

/opt/gridview//pbs/dispatcher/mom_priv/jobs/14923.admin.SC: line 6:
$logfile: ambiguous redirect

I need help on how to resolve this.
Thanks.

This most likely means your variable $logfile does not exist.

At a guess, the variable $logfile is not properly defined (i.e. it hasn't been defined at all). What is the output of the following command:
echo "logfile: ${logfile}"


text:
Fastq: how can I check if they are from DNA or RNAseq data?
I have (gave me) Illumina fastq files, which I want to use for variant calling, and I do not know if they are DNAseq or RNAseq data. How can I check this? I do not have any report or who to ask.
Many thanks

There are a few quick'n'dirty ways depending on the type of data.
In any case you want to align your files to a reference genome and then check the distribution of reads, either on a genome browser or with tools such as RSEQC
which calcualtes the fraction of reads aligning to exon, intron, intergenic etc.
RNA-seq, if you use a standard aligner such as bowtie2 or bwa, and you have a higher species that splices its RNA, then you should see most reads aligning to exons, and quite a fraction of reads that span a splice junction being unmapped (because the exon-intron-exon gap is large and non splice-aware tools do not "bridge" this without tremendous drop in mapping quality so the reads goes unmapped).
DNA-seq (if whole genome) should have a somewhat even coverage across the genome on the global scale. Targeted assays such as ATAC-seq should give distinct peaks that you can easily see in a genome browser just by eye.
If it is exome-seq (so targeted DNA-seq via exon pulldown) then it might get tricky because also most reads (like in RNA-seq) should align to exons, but you should not have spliced reads so the 5' part of the read starting in exon1 and 3' part ending in exon2. Instead many reads should also cover the exon-intron boundary.
As said, I would simply align it and check on a browser, I could imagine one can make a sophisticated guess by just looking at the alignments by eye.
A different topic is reliability, because I personally would not touch a dataset with unknown origin, because that also means you have no idea about the protocols that were used for library prep and what the sample actually is, so I wonder which question your analysis is going to answer.

Very nice post @ATpoint

Honestly, ask who made them.  If the reads are from an organism without splicing, it will not be trivial to figure it out with 100% certainty.
If you are calling variants, I don't think it matters much, except that higher organisms have such a small % of their genome transcribed, it should be obvious very quickly if you have the desired coverage across the entire genome, or just the transcripts.


text:
How to use combat in order to remove batch effects?
I have RNA seq data and I need to use combat to remove the batch effects. Somehow when I run it, it isnt actually doing anything.
The code:
dds <- DESeq(dds, betaPrior=TRUE)
normalized_counts <- counts(dds, normalized=TRUE)
log2 = log2(normalized_counts+1)

modcombat = model.matrix(~dex, metadata)# - with metadata being a variable containing treatment or control (under dex column) and batch and name of each patient.

com<-ComBat(log2, metadata$Batch, mod = modcombat)

Its supposed to be 4 different batches, but in the com variable I can see that the values have stayed the same as in log2.
What could be wrong? Would appreciate any help!

side note: combat is a linear batch correction method and does not perform optimally on all datasets

Are you completely sure that you can make the model matrix with a mix of column names as strings and object names like that?

No I am not sure about that, I havent found a good explanation for how to make the model matrix. I suspect my problem is there maybe.


text:
How to use CD-HIT to filter protein sequences dataset by certain similarity threshold?
I am looking at using CD-HIT to efficiently filter my protein sequences dataset by a similarity threshold of 70% (cut-off). More precisely, what I want to achieve is that for all the remaining sequences after filtering, all pairwise sequence similarity scores that can be computed are less than 70%.
Does that mean that I should run CD-HIT on my dataset with similarity threshold=70% (rest of the settings at default values), and then just keep only the "representative" sequences from the resulting cluster files? I have read the guide (http://weizhongli-lab.org/lab-wiki/doku.php?id=cd-hit-user-guide) but I am still not sure of how to use the output.
Linking to the above, how would CD-HIT handle this scenario below or similar ones like this:

assuming seq. similarity threshold=70%
sim.(seqA, seqB) = 80%
sim.(seqA,seqC) = 90%
sim.(seqA,seqD) = 65%
sim.(seqA, seqE) = 60%
sim.(seqD, seqE) = 50%
=> here, would the method keep just seqE, seqD in the output as sequences satisfying the threshold?

Thanks!

I don't follow your question. The sequences seqD and seqE do not fall into the seqA|seqB|seqC cluster, nor do the fall into the same cluster, so you will have three clusters, with the type sequences seqA,seqD and seqE. But that's not your question is it?


text:
add task to nextflow and resume
in test.nf (dsl2) I have
   ...
   workflow{
       a = Channel...
       processA(a)
       processB(processA.out.emitA)
       processC(processB.out.emitB)
    }

The pipeline completed successfully.
Now I want to add processC_2 that takes the output of processB.
How do I run the workflow to bypass processA execution?
I added processC_2 to test.nf and then did:
nextflow run test.nf -resume run-name

But it doesn't run the newly add process.
How do I achieve my goal?
Thank you!

It is actually very hard to answer this using only the information you have provided. It'd be necessary to actually see your test.nf file and read the processes. The execution seems right and I'd guess that the processB may be not emiting any output as shown in https://www.nextflow.io/docs/latest/process.html#outputs.
It is also possible that when using the -resume together with a specific run-name you may be recovering the state of that specific run as shown here https://www.nextflow.io/blog/2019/demystifying-nextflow-resume.html. You can try to execute it without the run-name or without the -resume flag so you can check whether the problem is with the emition and channelling of outputs between processes, or with the -resume flag.


text:
Align the HIV-1 Gag protein with the Gag protein of Visna virus
Align the HIV-1 Gag protein with the Gag protein of Visna virus
Background Visna virus is a lentivirus causing encephalitis in sheep
The problem is align the follow two protein Gag of viruses HIV-1 (GenBank KU749412.1) and Gag Visna virus (GenBank L06906.1) that are:
"MGARASILRGEKLDTWERIRLRPGGKKHYMLKHLVWASRELDRFALNPGLLETLEGCKQIIKQLQPALQTGTEELRSLFNTVATLYCVHAGIPVKDTKEALDKIEEEQNKCQQKAKQAEAAAAGKISLNYPIVQNLQGQMVHQAISPRTLNAWVKVIEEKAFSPEVIPMFTALSEGATPQDLNTMLNTVGGHQAAMQMLKDTINEEAAEWDRLHPVHAGPVAPGQMREPRGSDIAGTTSTLAEQIAWITGNPPVPVGEIYKRWIILGLNKIVRMYSPXSILDIRQGPKEPFRDYVDRFFKTLRAEQATQDVKNWMTDTLLVQNANPDCKTILRALGPGATLEEMMTACQGVGGPGHKARVLAEAMSQTTNAIMMQRSNFKGPKRTIKCFNCGKEGHLARNCRAPRKKGCWKCGKEGHQMKDCNERQANFLGKIWPSHKGRPGNFLQSRPEPTAPPAESFRFEETTPAPKQEQKDKEPLTSLKSLFGSDPLSQ"
"MAKQGSKEKKGYPELKEVIKATCKIRVGPGKETLTEGNCLWALKTIDFIFEDLKTEPWTITKMYTVWDRLKGLTPEETSKREFASLQATLACIMCSQMGMKPETVQAAKGIISMKEGLHENKEAKGEKVEQLYPNLEKHREVYPIVNLQAGGRSWKAVESVVFQQLQTVAMQHGLVSEDFERQLAYYATTWTSKDILEVLAMMPGNRAQKELIQGKLNEEAERWVRQNPPGPNVLTVDQIMGVGQTNQQASQANMDQARQICLQWVITALRSVRHMSHRPGNPMLVKQKNTESYEDFIARLLEAIDAEPVTDPIKTYLKVTLSYTNASTDCQKQMDRTLGTRVQQATVEEKMQACRDVGSEGFKMQLLAQALRPQGKAGQKGVNQKCYNCGKPGHLARQCRQGIICHHCGKRGHMQKDCRQKKQQGNNRRGPRVVPSAPPML"

I propose that the better alignment method is the one that maximize the function
f(string1, string2)=(number of letter equal between strign1 and string2)/(common length of strings).

With that definition I found f(s1,s2)=0.25760 (and 0.2949 between POL protein) for the two strings above.

What score might others find?
Is that definition one of biological significative?
Is there is a better one, what definition for align proteins?
Are Visna and HIV viruses relatives?

There are a number of scoring methods. Maybe review this page for an overview and compare different scoring metrics: https://en.wikipedia.org/wiki/Sequence_alignment

Both are lentiviruses and env is a fast evolving gene. What you are looking at here is the location of indels and this will vary between alignment algorithms.
There will always be an element of error in attempting to align over this genetic distance. It depends what you end goal is, if it is to maximise homology than okay, however optical enhancement of alignments are common if your downstream analysis is phylogenetics. If it is protein structure homology style modelling, where you have the protein structure of HIV env experimentally determined mulitple times you want to use it to help estimate that of the cat lentivirus then optimising a scoring metric is okay.
What you have to understand is there is no 'right' method, in phylogenetics this doesn't matter as much because indels are excluded in the subsequent calculation, so it is less sensative to indel positioning errors and optical enhancement is used. You can also simply exclude regions that appear 'garbled'. For protein homology modeling this is a big problem, because the indels will make drastic changes to the structure prediction. So the only approach here in my opinion is to produce a separate structural prediction for every variation in indel alignment you have obtained. In the latter context getting the alignment right is everything and I suspect thats the output here, because its a classic surface antigen to understand, which will readily crystalise.

yes, only it is Gag...  env protein would be the next... i want just check what are the origin of hiv


text:
Prokka error message
I am trying to run Prokka but the following error message appears:

Could not run command: makeblastdb -hash_index -dbtype prot -in /cluster/software/prokka/1.13.7-gompi-2019a/db/kingdom/Archaea/sprot -logfile /dev/null

When I removed the --setupdb, the following error message appeared:

The sequence databases have not been indexed. Please run 'prokka --setupdb' first

after that, I changed my script to
prokka --setupdb --kingdom Bacteria -o $OUTPUTFILE.txt o- $OUTPUTFILE.tsv o- $OUTPUTFILE.tbl o- $OUTPUTFILE.faa o- $OUTPUTFILE.gff /cluster/work/users/ahmedbar/191129_M01132.Project_Petersen-DNAlibs2-2019-11-06/prokka_A_C/$SAMPLE_ID_1

and still not working and the following error message appeared

Could not run command: makeblastdb -hash_index -dbtype prot -in /cluster/software/prokka/1.13.4-foss-2018b/db/kingdom/Archaea/sprot -logfile /dev/null

What is the output to which -a makeblastdb?

How did you install prokka?

It is installed by the admin. I work through the university server.

i think the answer to Ram RS's question will be crucial. what happens when you try to run makeblastdb -hash_index -dbtype prot -in /cluster/software/prokka/1.13.4-foss-2018b/db/kingdom/Archaea/sprot ? I wonder if you don't have write permissions to that directory, so it's failing. I would suggest contacting your admin, who will have a better idea of your setup than we have.


text:
Biopython SeqIO check input file
Hi I am trying to learnt python3 and Biopython,
I am trying to check imported fasta file before processing so far using:
sequences = SeqIO.parse("fake_fasta_file", "fasta")

if not any(sequa):
        print('\n  SOMETHING IS WRONG WITH YOUR INPUT FILE \n')
        exit()

works but as soon as I enter a '>' in my fake_fasta_file , as an example
nihvivicòvo-l
,mnjobùJN
, Békibpç
, mnkn

nbigp
ù
mlkn

mlnù
>  ftfyfiuuogpigèi

any(sequences) returns True even if:
for seq_record in sequences:
       print(seq_record)

doesnt print anything
Any better way to check for real fasta file in Biopython ?
Or maybe it my fault using any() on the iterator instead of creating
a list of seq_record and checking on that ?

any() loops over the iterator until the first truthy item is found or the end is reached. The iterator is not reset to the start after this. Using a list will change that behaviour since a list is not exhausted like a iterator. So, remove your any() call and you will see a sequence record printed. Your 'fake_fasta_file.fa' is actually recognised by Biopython as valid fasta:
>>> SeqIO.read('fake_fasta_file.fa', 'fasta')
SeqRecord(seq=Seq(''), id='ftfyfiuuogpigèi', name='ftfyfiuuogpigèi', description='  ftfyfiuuogpigèi', dbxrefs=[])

Biopython ignores lines before the first line beginning with '>' (as possible comment lines) and then starts collecting records, in this case an empty sequence record is still a record. This behaviour is deliberate from Biopython, although you could argue if it's best or not, it's unlikely to change for historical reasons.
What is real fasta? (philosophical pause) If you want to validate a particular definition of FASTA format then you can write your own parser depending on the exact specification and rules you want.
Edit, e.g.
from Bio import SeqIO

records = SeqIO.parse("myfile.fa", "fasta")
first = next(records, None)
if first is None:
    print("No fasta records")
elif first.seq == "":
    print("Empty sequence")
elif all(c not in "ATCGN" for c in first.seq.upper()):
    print("Non ATCGN characters")

etc.

thank You very much indeed for the very informative and detailed answer. Given I am really a newbie is there any chance you could point me to pre-made fasta validator I could borrow for my script ? (i was thinking about valid if >header than sequence on one or multiple lines and a way to check the characters used in the sequences)

@pippo1980 see edit (untested)

cheers !!!!!!!!!


text:
Repeat analysis on eukaryotic assemblies
I have a hundred insect genomes and I'm looking for repeated regions along these assemblies.
First of all, I thought of using RepeatModeler, only the most familiar among you know that it is not at all optimized software that requires a lot of resources and calculation time, which is therefore not feasible in my case.
One possible solution Repbase is not free.
I thought for example to use a blastn approach by using in query the assemblages and in database of nucleotide databases containing ETs, which would be an approach maybe less sensitive but much faster! Only I can't find a database in ET representative of eukaryotes
I have two questions:

Is there an alternative algorithm that is fast enough to process a large number of insect genomes?
Is there a database of representative ETs to enable a blastn style solution?

Thanks a lot.

when you say "ET" do you mean "TE" ("transposable element")? If not, what does that abbreviation stand for?

Yes Transposable element sorry


text:
Loop Pymol command get_area
I want to select multiple residues from different PDB files loaded on the same pymol session.
select resi LIG
get_area sele

How do I loop the above command to enable the command below:
select ligands,resn LIG

from https://www.blopig.com/blog/2019/08/how-to-iterate-in-pymol/
Example 1: Print the names of all the selected objects
You’ve selected several objects in PyMOL, and you want to see their names. By default, a new selection object called “(sele)” is created. We can create a Python variable called seleobjs, say, and then loop over them:
seleobjs = cmd.get_object_list('(sele)')
    for obj in seleobjs: print obj

Example 2: Disable all the selected objects
If you still have your selection, “(sele)”, you can disable them (effectively turn them off on the right-hand sidebar):
seleobjs = cmd.get_object_list('(sele)')
    for obj in seleobjs: cmd.disable(obj)

Example 3: Delete all the selected objects
You’re done with the objects in your selection, “(sele)”: you might be able to delete them using “delete ”, but if you can’t write a wildcard that covers all the names of the objects, you can iterate over a selection and delete them:
seleobjs = cmd.get_object_list('(sele)')
    for obj in seleobjs: cmd.delete(obj)

from https://pymolwiki.org/index.php/Get_Area
get_area calculates the surface area in square Angstroms of the selection given. Note that the accessibility is assessed in the context of the object(s) that the selection is part of. So, to get the surface areas of e.g. a component of a complex, you should make a new object containing a copy of just that component and calculate its area.
please let us know if it worked

Thanks, this code is working, but I am still not able to calculate the area of selected residues. My main hurdle is how to use the get_area pymol command and incorporate in the loop.

@shome edited my answer

I typed :
seleobjs = cmd.get_object_list('(sele)')
for obj in seleobjs: get_area obj
It is showing me invalid syntax..

cmd.get_area(str selection="(all)", int state=1, int load_b=0)     try with Parentheses (Round Berackets)

showing invalid syntax

both get_area(obj) and cmd.get_area(obj) ?

Ok went back to this question having pymol within a hand's reach this time. Question is not really clear to me so I am posting a reasonable (to me) version of the question into a new code snippet;
I am using two PDBs containing ATP as ligand '3dgl.pdb', '2p09.pdb'
code:
import pymol

list_pdbs = ['3dgl.pdb', '2p09.pdb']

pymol.finish_launching()

for _i in list_pdbs :

    pymol.cmd.load( _i, _i.rstrip('.pdb'))

objs = pymol.cmd.get_object_list(selection='(all)')

print('structure loaded : ' , objs)

pymol.cmd.remove('resn HOH')

for _i in objs:

    pymol.cmd.create(_i+'_ATP_copy' ,  _i + ' and resn ATP')

    # print('after create ')

    area = pymol.cmd.get_area(_i+'_ATP_copy', quiet = 0)

    pymol.cmd.show('dots' , _i+'_ATP_copy')

    print(_i+'ATP_copy' + ' area --->  ' + str(area))

output:
structure loaded :  ['3dgl', '2p09']
 cmd.get_area: 1013.704 Angstroms^2.
3dglATP_copy area --->  1013.7040405273438
 cmd.get_area: 1013.704 Angstroms^2.
2p09ATP_copy area --->  1013.7040405273438

picture of the ATP surfaces in Pymol itself,

you can play with  cmd.get_area(...) parameters see https://pymolwiki.org/index.php/Get_area


text:
How can I get properties from each one letter amino acid?
Im using python and wish to get properties from each one letter amino acid?
I mean
M(Met) => Density:1.340, MW:149.21...
The RDkit methods atom.GetSymbol(), atom.GetDegree() permit this and I am looking for an equivalent solution in Python.
Are there any suitable package, for example in Biopython?

Check ProteinAnalysis from Bio.Seq.Utils.ProtParam in Biopython.

RDKit is a Python solution.


text:
How to extract clusters from a phylogenetic tree?
I am using Biopython Phylo and RDkit Cluster to obtain a UPGMA tree from a distance matrix of 5k x 5k entries. How can I extract the taxa names within a given clade(s)?
For example, Bio.Phylo can return plain-text representation of a tree, so I think I can parse it by the number of tabs:
Phylogeny()
     Clade()
         Clade(branch_length=0.06)
             Clade(branch_length=0.102, name='A')
             Clade(branch_length=0.23, name='B')
         Clade(branch_length=0.4, name='C')
...

or, ASCII format:
          _____________ A
  _______|
_|       |_______________________________ B
 |
 |_______________________________________________________ C

Which, I think is more difficult to parse but possible.
Either way, those are not established methods for this kind of task, so I would like a hint.

There are two packages that do this well, ETE3 and dendropy. I will put a ETE3 solution forward.

The example given a 3 taxa tree, isn't an ideal example. What is required is a raw treefile in a common format (nexus, philip, newark). The example below uses ETE3 for a 5 taxa tree A,B,C,D,E. Its ROOTED - very important to know the direction.
ETE3 is providing a means of traversing down the tree from the root, however the coder needs to interrogate the path it has taken.

Tree file imported
Ladderizing ... so you know the horizontal direction of the tree (see *)
Ask the ETE3 to search down the nodes ... each child is a node from the root.
Parse the output and identify the clade.

The two outputs from point 3 are either it encounters a node - in this case it returns all the children below it (taxa names). The alternative is that it encounters a terminal node or a leaf, like a species name or isolate. In this case the taxa name is returned.
The caveat is the terminal node (taxa label or leaf) is always bifurcation and the taxa names are split into [1] and [0] on the ETE3 object (see code) ... you have to interrogate both.
Thats what the OP is looking for and this part of the code is marked 'Part 1'
Parsing Part 2 is more complicated, because usually the coder wants the clade a particular taxa (or leaf) belongs to. Essentially, the coder needs to ask what the leaf name is, if the target leaf name is found then the clade is then identified, i.e. the other stuff next to it.
The code is below. I agree if you don't think trees, this is a rushed explanation.
from ete3 import Tree

# Part 1
tree = Tree(intree, format=1)
tree.ladderize(direction=1)
clade1_ = tree.children[0].children[0].children[0].name
leaf1_ = tree.children[0].children[0].children[0].children[0].name
leaf2_ = tree.children[0].children[0].children[0].children[1].name

# Part 2
if leaf1 == 'A' or leaf2 == 'A':
    CladeA = clade1_
if leaf2 == 'B' or leaf2 == 'B':
    CladeA = clade1_

*, ladderizing that can be important for classifying [0] and [1] in ETE3 ... its a weird one in trees 'cause it shouldn't matter


text:
Why do I obtain different output results with blast vs awk commands
I have an awk command that identifies 30 pb from two multifasta files.
When I used two input files:

E.g. 100 sequences each, I get the same result with the awk command and blastn command.
E.g. with 1,000 sequences each I get the same result with both commands.

Question However, when I increase the number of sequences to 10,000 I get different outputs, any idea why this could happend?
Command-#sequences-matches
awk      -100-       1
blastn   -100-       1
awk     -1000-      73
blastn  -1000-      73
awk     -10000-     8384
blastn  -10000-     8280

awk command:

BEGIN { wid = 30 }

sub(/^>/,"") { hdr=$1; next }

NR == FNR { a[hdr]=$0; next }

{

    for ( hdrA in a ) {

        strA  = a[hdrA]

        lgthA = length(strA)

        for ( idxA=1; idxA<=(lgthA - wid); idxA++ ) {

            substrA = substr(strA,idxA,wid)

            if ( index($0, substrA) ) {

                printf "[%s,%s]\n", hdrA, hdr

                break
            }
        }
    }
}

awk -f tst.awk file_with_100seqs file_with_100seqsv2

blastn command:

blastn -query file_with_100seqs -subject file_with_100seqsv2 -out blastn_output_seqs1vs2_100.txt -word_size 30 -perc_identity 100 -outfmt "10 qseqid sseqid"
The input files look like this:
#File1
>1
TGATTGCATAACCACTTAACATCTTGTTTTATCTAAATAAAATTAAGCATGTTATCTTTTTGGGGCACTCCTGGGGCAGTAGATGCCAGTTGTTGATTCAGTATATCTACTTGTGACTGGTTATTATCCCGATTTTTTTAGTTTTAAGGTGTTGACATAGCCATCCATGCTCCATATACTGTATAGACCATCTGAGCGTT
>2
TGGGAAAACAGCATTCAGCGGTGGCTTATTCCTGCTAAGGATGTTGGCCGCATTCATGCTGAGCACAACCTCGACGGCCTGCTGAGGGGCGATTCGGCATCCCGCGCTGCCTTTATGAAGGCAATGGGAGAGGCAGGGCTACGCACCATCAACGAGATGCGACGAACGGACAACCTCCCGCCATTGCCGGGTGGCGATGT
>3
GAAATGGGAACCGCGAACATGCCTGCACATCCGTTTGTGCGACCCGCTTACGATACTCGCGAGGAAGAGGCCGCCAGCGTCGCCATTGCCAGGATGAATCAGGCTATTGATGAGGTATTGAGCAAGTGAATGAAGATAATATCTACGCCTTGCTTTCTCCCCTGGCAGAAGGACGGGTATATCCCTATGTTGCGCCATTA

#File2
>1
TGATTGCATAACCACTTAACATCTTGTTTTATCTAAATAAAATTAAGCATGTTATCTTTTTGGGGCACTCCTGGGGCAGTAGATGCCAGTTGTTGATTCAGTATATCTACTTGTGACTGGTTATTATCCCGATTTTTTTAGTTTTAAGGTGTTGACATAGCCATCCATGCGGGAAGGTGCAGCATAATGTGCTTTGGATT
>2
TGAGTGCCCCATTTGTGAAGCAATAAAGTTCGGGTTCGCGCCAGCGGCAAGCGCCCAGCATGCACCGATTTTTTTAGTTTTAAGGTGTTGACATTAGGTATGTCGGGACTGGTATGCTTTCCTGTGTCGCAGCCCGGCGCGTCTCAATGCAGATTCCCATATCCTGTTCATCCATATACTGTATAGACCATCTGAGCGTT
>3
TACCTGAGCGATCGGTAATTTGCGGATTGAAGACAAAGGTGCAGGAATGAGTTTTTGTACGACCGTATTCGCGCAGCTTTACTTCAATTTTGTGCTGTTTGCTCAGCTTCGTGAAAGAGGCCTGACTTTTTAAAGCATCAATTGCTGGCTGCACAAGATGTATCACCCTGTCGGTTCCTGCCTGGGTTTTCGGCAGGGTG

Can you please tell us what you are trying to accomplish? It is not clear what your input data are, what you are searching for, and what the intended output is.


text:
What is the origin of HIV1?
What is the origin of HIV1?
This is the alignment of a small part of the POL protein of the HIV1 virus GenBank KU749412.1, with the POL protein of the Visna virus GenBank L06906.1
IKKKDSTKWRKLVDFRELNKRTQDFWEVQLGIPHPAGLKKKKSVTVLDVGDAYFSVPLDESFRKYTAFTI--PSTNNETPGIRYQYNVLPQGWKGSPAIFQSSMTKILEPFRAQNPEIVIYQYMDDLYVGSDLEIGQHRAKIEELREHLLRWGFTTPDKKHQKEPPF--------LWMGYELHPDKWTVQPIQLPE--KDSWTVNDIQKLV
|||| | ||| | ||||||| | |  | ||| ||| ||  || || || |||||  || |  | || |    |  ||  |  ||   ||||||| |||  |  | |||       | |    |||| | ||||  |      || |                                | |||| ||  |   |||      | |  ||||
IKKK-SGKWRMLIDFRELNKQTEDLAEAQLGLPHPGGLQRKKHVTILDIGDAYFTIPLYEPYRQYTCFTMLSP--NNLGPCVRYYWKVLPQGWKLSPAVYQFTMQKILRGWIEEHPMIQFGIYMDDIYIGSDL--GL-----EEHRGIVNELASYIAQYGFMLPEDKRQEGYPAKWL-GFELHPEKWKFQKHTLPEITEGPITLNKLQKLV

Are these matches,

By chance or convergence,
HIV1 and Visna virus share a recent common ancestor or,
the common ancestor was closer Visna virus from where HIV descended?

Could it be that HIV1 really comes from some kind of experiment involving Visna virus, as I read on the Internet?
The disease caused by the Visna virus is (100% deadly see https://pubmed.ncbi.nlm.nih.gov/16425963/ ). Could there be parallels to sars2?
But look at it: the Visna virus is transmitted by air the first symptoms are like pneumonia (or a brain infection) and otherwise are chronic progressing slowly like SARS-CoV-2 infection resulting in long Covid.
Symptoms are also felt after many days that the virus has come into action, this seems to me not to happen for coronaviruses or for flu viruses but for HIV and Visna virus.
Related links:

https://pubmed.ncbi.nlm.nih.gov/8289690/

https://pubmed.ncbi.nlm.nih.gov/16425963/

(preprint sars2 as retrovirus) https://www.biorxiv.org/content/10.1101/2020.12.12.422516v1

https://www.gao.gov/assets/100/91319.html

one 1966 HIV1 virus in https://www.pnas.org/content/117/22/12222  MN082768

They are both lentiviruses and share a distant common ancestor.
HIV-1 and HIV-2 are descendents from simian immunodeficiency virus (SIV). The following tree shows the relationships very clearly, from Wertheim and Worobey (2009) Dating the Age of the SIV Lineages That Gave Rise to HIV-1 and HIV-2 PLoS Comp Biol. here. The evolution of HIV-1 is heavily mixed with chimanpzee (cpz) SIV, and there are regions of highly anomalous homology juxtaposed with regions of notably lower homology.

and where it came from "SIVcol" and what is his genome?

Colobus monkey https://retrovirology.biomedcentral.com/articles/10.1186/1742-4690-10-107

so you say visna->siv->HIV1 & HIV2?  what about the chain is reversed visna->HIV1 & HIV2 ->SIV

No, he's saying OriginalLentivirus->SIV/HIV and OriginalLentivirus->visna.

So seems to me understand genome SIVcol can be find in GenBank with number
AF301156.1 (SIVcolCGU1), with other SIVcol gag protein KF214242.1,
with other SIVcol partial genome KF214241.1, the one base to all SIV Tree.
SIVcol the one all other SIV they claim came from, right?

@mgkrebbs no in the link they claim siv->HIV1 & HIV2  so first was siv virus than HIV1 and HIV2 they say

@RosLuP, yes that's true too, and well known, but the point is there is no evidence that HIV-1 is anymore related to visna than to other lentiviruses, such as the lentiviruses of cattle and the cat family.


text:
Is there a straightforward way to get mismatches/indels in a BAM file using pysam?
Ideally, I'd like to do something like the following:
f = pysam.AlignmentFile(filename, index_filename=index_filename)
read = list(f.fetch('Chromosome', 1, 10000))[0]
read.get_variations()

And have that output something like:
[
  (1, 'A', 'C'), # An A -> C mismatch at position 1
  (15, 'I', 'CGGG'), # An insertion of the string 'CGGG' starting at position 15 of the read
  (32, 'D', 4), # A deletion of 4 base pairs starting at position 32
]

I know that it's possible to figure out all of those values by parsing the CIGAR string and MD tags but I'm wondering if there's a more intuitive interface.

You're looking for pileup, which is the htslib (and thus samtools/bcftools) method for finding variants. See this section of the pysam documentation. You may want to peruse about the output as well. The gist is:
import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb" )
for pileupcolumn in samfile.pileup("chr1", 100, 120):
    for pileupread in pileupcolumn.pileups:
        if not pileupread.is_del and not pileupread.is_refskip:
            # query position is None if is_del or is_refskip is set.
            ...handle mismatches...

Thanks for the response! I think this would work but it's less than ideal since I would be iterating over each read multiple times, no? I'm wondering if it's possible to iterate over the list of reads that span a region and retrieve that same information without going over pileup columns.

Also, this bit of code for pileupread in pileupcolumn.pileups crashes the Python kernel when running inside of an IPython notebook :-/


text:
Merging plink files with divergent allele codings
I have a number of plink files from different arrays (mixed Illumina and Affymetrix) that need to be merged. However, these use different allele codings, so trying to merge without fixing this causes an error (3+ alleles). One can fix this by using (on each file) --recode 12 to a PED file, and then converting back to BED (--make-bed) and then merging (--merge-list). Doing this, I was able to merge the files without explicit errors. However, this method is suboptimal, I think, because a minor allele in dataset 1 might not be the same one as in dataset 2, since this depends on the people in the dataset.
For instance, rs746969. Illumina 1MDuo:
10  rs746969    1.00698e+06 100008620   1   2

Affymetrix 6.0:
10  rs746969    0   100008620   A   G

Alternatively, I could keep the 12-encoded merged data, and attempt to convert back to ACGT. Which is the best approach?

Yes I would advise against doing the first option, since like you say, you don't know the minor allele is going to be the same between the different datasets. When merging two different datasets, it is very important that you make sure there aren't any strand flips and that the alleles line up. I would reccomend making sure that both datsets are in the ACGT format and then using something like genotype harmonizer to check for strand flips.


text:
How to calculate protein abundance from an .mzTab file?
I have an .mzTab file, with PRT and PSM values.  Does anyone know of any specific functions that I could use to calculate protein abundance values from this data?
More details
This is a summary of what the mzTab file is.
The data I have has two main sections within the file - a PRT section and a PSM section.
The PRT section has protein information.  It has the following column headers:
 [1] "accession"                        "description"
 [3] "taxid"                            "species"
 [5] "database"                         "database_version"
 [7] "search_engine"                    "best_search_engine_score[1]"
 [9] "search_engine_score[1]_ms_run[1]" "num_psms_ms_run[1]"
[11] "modifications"

Where num_psms_ms_run is the count of the total significant PSMs that can be mapped to the reported protein.
The PSM section has peptide spectrum match information.  It has the following column headers:
 [1] "sequence"
 [2] "PSM_ID"
 [3] "accession"
 [4] "unique"
 [5] "database"
 [6] "database_version"
 [7] "search_engine"
 [8] "search_engine_score[1]"
 [9] "search_engine_score[2]"
[10] "search_engine_score[3]"
[11] "modifications"
[12] "retention_time"
[13] "charge"
[14] "exp_mass_to_charge"
[15] "calc_mass_to_charge"
[16] "spectra_ref"
[17] "pre"
[18] "post"
[19] "start"
[20] "end"
[21] "opt_global_mzidentml_original_ID"
[22] "opt_global_cv_MS:1002217_decoy_peptide"
[23] "opt_global_cv_PRIDE:0000091_rank"

More information about the column definitions can be found here.

Hi @icedcoffee please do supply more detail for example I'm not sure what the acronyms are, I'm pretty sure your question is fairly easy whether directly or indirectly by customised code.

Thank you, I've updated the question to include more information


text:
Which one is a more convenient assembly?
I have developed a software for de novo genome assembly. Its performance varies gradually according to how much data you employ. At initial stages it often produces contigs that look like that when aligned to a reference genome:

Initial-stage s-aligner assembly of SRR8357422
Genome fraction:98,8% -
Largest Alignment: 61.544 -
NGA50: 59.409 -
LGA50: 2
If I add more data, the contigs end up being larger, but that obviously takes more time. My question is, how does this result compare to, for example, one like that?

SPAdes assembly of SRR8357422
Genome fraction:93,46% -
Largest Alignment: 128.391 -
NGA50: 128.391 -
LGA50: 1
Note how you can construct a contig covering almost the full reference genome assembling green contigs in the first case. Therefore, a blast search of genes will likely let you find more genes given that it covers more genome fraction.
Apart from determining what option is more convenient, I would also like to know how much more convenient is to have both results. Or in which cases it's more convenient.
Last (but not least!), is there a way to generate automatically larger contigs from these overlapping contigs in assembly 1? I have tried using miniasm with no success. Even minimap2 seems not to find overlaps well enough. I am developing my own algorithm also for that but would like to know if there is another way.
[EDITED on 17/12 to add another example with real data for both assemblies]

Hi, could you please clarify your question a bit more; it's unclear to me. Asking "which one is a more convenient assembly?" suggests to me that there should be at least two comparisons, but you've only provided a figure for one. There is also additional information needed about the assembly and reference genome (i.e. are you doing kmer-based assembly; is the genome highly-repetitive; does it have Ns in the sequence).

Hi @gringer, it's not a question for a single case. It's something that happens in different assemblies in different circumstances. It's a generic question. The example in the image is real and corresponds to a run in PRJEB32127 but I have not assembled it with another software yet. I mean, in case I obtained an assemblie with a single contig but lower coverage (for example a contig 100k long) , which one of the two assemblies would be more useful.

Maybe it's more clear with a more extreme case. Suppose it's a bacterial assembly (~2Mbp) and SPAdes get a contig with 3Mbp bit just covering 90% of the reference (the rest is chimeric) while the otrher assembler gets lots of contigs around 200k bp long but covering 99%. I'd like to know how important is to get larger contigs or if coverage and percentage of genes covered is more important.

I will lsee if I find a real case with real data...

I guess it's "genome fraction" where I said "coverage"...

If it's a general question about assembly algorithms, it won't really fit well into the StackExchange Question/Answer format. Bioinformatics Stack Exchange works best for single, specific issues. I recommend posting your question on the bioinformatics subreddit.

Yes, makes sense. I will post it there. I will see if I find an example clear enough too and edit the answer or write a new one.

I asked that on Reddit. I didn't get an answer to all my questions but I got an answer for how to assemble the overlapping contigs.
If I use Flye I can assemble these and get larger contigs, better NG50, and higher genome fraction.

s-aligner + Flye assembly of SRR8357422 Genome fraction:95,95% - Largest Alignment: 134.103 - NGA50: 134.103 - LGA50: 1
Given that the data in the first assembly apart from being more useful itself in my opinion, is objectively more suitable to generate a better result in all senses with a simple additional step that takes a few seconds, my only possible conclusion is that it is overall a better result and therefore this sort of result will often be preferable.
Of course, there are other aspects to consider, like error rate, that are not considered here as it wasn't part of the question (mainly because I think Quast is not adequate for calculating it). But BTW, mismatches per 100kbp are also lower.
[EDIT] I want to add that this is not granted that in all cases Flye will generate a clearly better assembly.


text:
FastQC and Trimmomatic in Galaxy?
I am independently working on data retrieved SRA database, paired-end data as separate inputs.
The proceedure I followed is,

After running FastQC using Galaxy, the majority of the modules have failed.
I tried Trimmomatic with default settings (Avg quality =20, number of  bases to avg across = 4), it resulted in R1 and R2 paired and unpaired (4 total outputs).
I ran FastQC again, which rendered some changes (per base sequence quality passed) but resulted in other modules (Sequence length distribution and per tile quality scores) getting a warning sign or failing.

My question is,

Is  there specific setting other than default that provides better results? What output should i use from Trimmomatic for fastqc again?

Hey, could you please be more specific about your data? What type of sequencing technology has been used? And what exactly are you trying to achieve? Are you trying to make all FastQCs modules pass the test ?

Also note that, depending on your data, you might not WANT all fastqc checks to pass. If you're using 16s amplicons of a limited set of different microbes you will always get fastqc warnings on some checks due to the nature of the data.

Thank you! its Illumina MiSeq paired end sequencing. I figured not all QCs can pass in this dataset. I am still confused to what output I have to move forward with the alignment?


text:
How to merge back the processes by sample after splitting by chromosome in Nextflow?
I have a pipeline that has an input list of samples. I would like to split those samples by chromosomes, process each chromosome apart, and then merge them back by sample name.
Here's an example:
samplesList = Channel.from(["sample1", "sample2"])

// A process that splits the samples by chromosome as an example
// In my actual case, I have a VCF that is split by `snpSift Split` command
// and generates files like [sample]_[chr].vcf
// In the example below, I only gave sample1
// But take into consideration that I have multiple samples

process bar {
   echo true

  output:
    file("*.tsv") into bar_output

  script:
  """
  echo 1 > sample1_chr1.tsv
  echo 2 > sample1_chr2.tsv
  echo 3 > sample1_chr3.tsv
  echo 4 > sample1_chr4.tsv
  echo 5 > sample1_chr5.tsv
  """
}

// Process each chromosome apart for each sample
// This means the output of the process foo
// is only ONE chromosome per sample
process foo {
  input:
    each x from bar_output  // to treat each file separately from the channel
  output:
    file("*.tsv") into xasxa
  script:
    """
    echo \\\$(cat $x)\\\$(cat $x) > sample1_chr\\\$(cat $x)_2.tsv
    """
}

// How should I write the next process to group back them together by sample name?
// I do not wish to do a collect() on channel foo because it will force
// all the processes of foo to finish first.
// In this example, foo would emit 5 different channels of chromosome for sample1
// This process should merge sample1 only when all chromosomes are processed
// independently of other samples.
process mergeChromosomesBySamples {
   input:
     ???
}

As the example above, the process mergeChromosomesBySamples should return something similar/close to
["sample1", ["/path/to/workdir/sample1_chr1.tsv",
             "/path/to/workdir/sample1_chr2.tsv",
             "/path/to/workdir/sample1_chr3.tsv",
             "/path/to/workdir/sample1_chr4.tsv",
             "/path/to/workdir/sample1_chr5.tsv"]
]

EDIT1: more explication for the process mergeChromosomesBySamples is given in the comments and as follows ...
The channel xasxa processes each chromosome apart, this means that the output of the channel xasxa contains only one chromosome of the sample like "/path/to/workdir/sample1_chr1.tsv" only. My goal is then to merge the processes out xasxa for each sample together.
I do not wish to use collect() in this case because it would have to force all the processes of foo to finish before. I want them to work independently for each sample.
Where each sample has its own VCF generated and I process them independently.
EDIT2: Here's a graph of what I am trying to achieve:

If I use .collect(), then all the processing on the chromosomes of all the samples has to finish first before being able to merge them, which is something I do not want.
I appreciate your help.

I was able to find a solution that answers my problem finally.
So, here are two different methods below, one which does exactly what I want while the other has to wait for all the previous processes to finish:
First method: waiting for previous process to finish using .groupTuple() like what @Pallie mentionned. This method is simpler than the second method. The problem is that it has to wait for the previous processes to finish. I want my samples to be treated in parallel. If I missed something with the logic of the method, please feel free to edit.
#!/usr/bin/env nextflow

// my samples
samplesList = Channel.from(["sample1", "sample2"])

// the process that is supposed to split the sample by chromosomes
process foo {
      echo true
      cpus 1
      memory '1 MB'
      executor 'local'

      input:
         val(x) from samplesList
      output:
         file("*.vcf") into foo_out
         val(x) into output_sample
      script:
      """
         echo 1 > ${x}_chr1.vcf
         echo 2 > ${x}_chr2.vcf
      """
}

//Combine the samples by chromosome and filter out only the samples that match the files output
foo_out.flatten().combine(output_sample).filter{ it[0].toString().contains(it[1].toString()) }.set{maybework}

// the process that is supposed to treat each chromosome apart
process bar {
      echo true
      cpus 1
      memory '1 MB'
      executor 'local'

      input:
         set file(y), val(x) from maybework
      output:
         set val(x), file("*chr*vcf") into bar_out
      script:
      """
         echo ${x}
         filename=\\\$(basename -- ${y})
         chromosome=\\\$(echo "\\\${filename%.*}" | cut -f2 -d'_')
         cat ${y} > ${x}_"\\\$chromosome"_"\\\$chromosome".vcf
         if [ ${x} == "sample2" ]; then sleep 20; fi
      """
}

// the process that merges back the chromosomes by sample
process mergefoobar {
      echo true
      cpus 1
      memory '1 MB'
      executor 'local'

      input:
      set val(sample), file(x) from bar_out.groupTuple()

      script:
      """
         echo $sample
         echo ${x}
      """
}

will output:

As you can see above, the process mergefoobar is still waiting for all the processes of bar to finish.
Second method: does not wait for the previous processes to finish, uses .join() for each chromosome... Yes, the downside is that you would have to write the join() for all the chromosomes. Additionally, if you use snpSift Split, one sample may generate 3 files out of 5 (because only 3 chromosomes had variants listed in them), you would need to create the other empty VCF files because the output excepts them.
#!/usr/bin/env nextflow

// my samples
samplesList = Channel.from(["sample1", "sample2"])

// the process that is supposed to split the sample by chromosomes
process foo {
      echo true
      cpus 1
      memory '1 MB'
      executor 'local'

      input:
         val(x) from samplesList
      output:
         file("*.vcf") into foo_out
         val(x) into output_sample
      script:
      """
         echo 1 > ${x}_chr1.vcf
     echo 2 > ${x}_chr2.vcf
         echo X > ${x}_chrX.vcf
     echo Y > ${x}_chrY.vcf
      """
}

//Combine the samples by chromosome and filter out only the samples that match the files output
foo_out.flatten().combine(output_sample).filter{ it[0].toString().contains(it[1].toString()) }.set{maybework}

// the process that is supposed to treat each chromosome apart
process bar {
      echo true
      cpus 1
      memory '1 MB'
      executor 'local'

      input:
         set file(y), val(x) from maybework
      output:
         set val(x), file("*chr1*vcf") optional true into bar_out1
         set val(x), file("*chr2*vcf") optional true into bar_out2
         set val(x), file("*chrX*vcf") optional true into bar_outX
         set val(x), file("*chrY*vcf") optional true into bar_outY
      script:
      """
         echo ${x}
     filename=\\\$(basename -- ${y})
         chromosome=\\\$(echo "\\\${filename%.*}" | cut -f2 -d'_')
         cat ${y} > ${x}_"\\\$chromosome"_"\\\$chromosome".vcf
         if [ ${x} == "sample2" ]; then sleep 20; fi
      """
}

// the process that merges back the chromosomes by sample
process mergefoobar {
      echo true
      cpus 1
      memory '1 MB'
      executor 'local'

      input:
      set val(sample), file(x), file(y), file(i), file(t) from bar_out1.join(bar_out2).join(bar_outX).join(bar_outY)

      script:
      """
         echo $sample
     echo ${x}
         echo ${y}
     echo ${i}
         echo ${t}
      """
}

Above, you see that the process mergefoobar has launched even if the previous bar has not finished yet.
I hope this also helps anyone in the future who might have this issue.

Use the groupBy operator. If the content of channel xasxa is:
["/path/to/workdir/sample1_chr1.tsv", "/path/to/workdir/sample1_chr2.tsv", "/path/to/workdir/sample2_chr1.tsv", "/path/to/workdir/sample2_chr2.tsv
You could probably do:
xasxa.groupBy{ it -> it.fileName.toString().split("_")[0] }
The result would be a channel with:
[sample1:[/path/to/workdir/sample1_chr1.tsv, /path/to/workdir/sample1_chr2.tsv, sample2:[/path/to/workdir/sample2_chr1.tsv, /path/to/workdir/sample2_chr1.tsv]]
This breaks if you have an underscore in your samplename before _chr*. More robust would be to initialise all channels with a sample value, so that each channel is in the form of [samplename, chromosome]. Then all you need is a groupTuple on your channel.

The channel xasxa processes each chromosome apart, this means that the output of the channel xasxa contains only one chromosome of the sample like "/path/to/workdir/sample1_chr1.tsv" only. My goal is then to merge the processes out xasxa for each sample together. Sorry if that wasn't clear in my explanation.

Still not exactly sure what you need, does a .collect() on  xasxa do what you wish?

I do not wish to use collect() in this case because it would have to force all the processes of foo to finish before. I want them to work independently for each sample.

Do you start with a single vcf that contains multiple samples?

No, each sample has his own VCF generated and I process them independantly. BTW, thanks for the proposed solution.


text:
How can I convert Protein ID to PDB ID?
I have a list of Protein ID such as Q9S9U6, A2XLL2, Q06402.
What name is this ID? UniProtKB ID?
And how can I get PDB ID from this list?
Is its PDB ID only one?
I'm using Python.
Get a list of PDB ids from a list of sequence protein number?
That method of this link gives me some IDs but not only one and usually empty.

I just wanna get an adjacency matrix from those IDs(Q9S9U6, A2XLL2, Q06402).

https://www.ebi.ac.uk/pdbe/docs/sifts/quick.html has the file you need for Uniprot <-> PDB

Not all proteins have knows structures, so it is expected that PDB IDs will not exist for all input UniProt IDs. Plus, some protein structures exist in multiple conformations, so there is no 1:1 mapping either.

Get their sequences and then search the pdb (by his API) for    the obtained sequences
https://search.rcsb.org/#search-services
...
sequence
This service employs the MMseqs2 software and performs fast sequence matching searches (BLAST-like) based on a user-provided FASTA sequence (with E-value or % Identity cutoffs). Following search targets are available:
pdb_protein_sequence: all current protein sequences in PDB...

I think this should help you : https://biit.cs.ut.ee/gprofiler/convert


text:
the model matrix is not full rank : This is a classic question which a biologist face without clear understanding of the model design
Saw this answer biostar. Tried to make my metadata as such but still the
"Error in checkFullRank(modelMatrix) : "

This is my coldata
dput(coldata)
structure(list(Group = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L,
6L, 6L, 6L, 6L), .Label = c("HSC", "Blast", "CMP", "GMP", "LSC",
"Mono"), class = "factor"), Mutation = structure(c(10L, 4L, 3L,
5L, 2L, 3L, 9L, 11L, 1L, 1L, 7L, 7L, 10L, 7L, 8L, 2L, 10L, 10L,
10L, 10L, 2L, 2L, 7L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 2L, 2L, 5L, 2L,
9L, 1L, 1L, 8L, 6L, 6L, 6L, 6L, 6L, 6L), .Label = c("DNMT3A-R882H",
"FLT3-ITD", "IDH2-R140Q", "KRAS-G13D", "MAU2-Q98H", "No", "NPM1-L287ins(TCTG)",
"NRAS-G13D", "STAG2-R614", "TET2-E1357", "TET2-R550"), class = "factor"),
    Satus = structure(c(2L, 5L, 4L, 6L, 3L, 4L, 9L, 10L, 1L,
    1L, 7L, 7L, 2L, 7L, 8L, 3L, 2L, 2L, 2L, 2L, 3L, 3L, 7L, 16L,
    16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L,
    16L, 16L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 12L, 12L, 13L,
    12L, 9L, 11L, 11L, 14L, 15L, 15L, 15L, 15L, 15L, 15L), .Label = c("Blast-DNMT3A-R882H",
    "Blast-E1357", "Blast-FLT-ITD", "Blast-IDH2-R140Q", "Blast-KRAS-G13D",
    "Blast-MAU2-Q98H", "Blast-NPM1-L287ins(TCTG)", "Blast-NRAS-G13D",
    "Blast-STAG2-R614", "Blast-TET2-R550", "LSC-FLT-DNMT3A-R882H",
    "LSC-FLT-ITD", "LSC-MAU2-Q98H", "LSC-NRAS-G13D", "Mature-Normal",
    "Progenitor-Normal", "Stem-Normal"), class = "factor")), row.names = c("Blast1",
"Blast10", "Blast11", "Blast12", "Blast13", "Blast14", "Blast15",
"Blast16", "Blast17", "Blast18", "Blast19", "Blast20", "Blast2",
"Blast21", "Blast22", "Blast23", "Blast3", "Blast4", "Blast5",
"Blast6", "Blast7", "Blast8", "Blast9", "CMP1", "CMP2", "CMP3",
"CMP4", "CMP5", "CMP6", "CMP7", "CMP8", "GMP1", "GMP2", "GMP3",
"GMP4", "GMP5", "GMP6", "GMP7", "HSC1", "HSC2", "HSC3", "HSC4",
"HSC5", "HSC6", "HSC7", "LSC1", "LSC2", "LSC3", "LSC4", "LSC5",
"LSC6", "LSC7", "LSC8", "Mono1", "Mono2", "Mono3", "Mono4", "Mono5",
"Mono6"), class = "data.frame")

dds <- DESeqDataSetFromMatrix(countData=df2, colData=coldata, design= ~ Group + Mutation)

This is the design I'm trying to incorporate. What i get is my design is nested in one of the column by looking at numerous example.
What i would like to do is I'm trying to run a LRT test on this dataset which works fine if i give "Group" as my design and then i run this
dds_lrt <- DESeq(dds, test="LRT", reduced = ~ 1)

which gives me diffrenitaly expressed across multiple groups. Now there are two question

Is it possible to include interaction such as in my case Group+Mutation and do the lrt test?
Or the above is conceptually wrong to do ?

what is the hypothesis you are trying to test?

I would like to know how these specific mutation affect gene expression across the cell types with respect to reference which is my "HSC" as its a stem cell . But from @ATpoint answer it seems I need to subset the data more over what i tired above perhaps is a wrong approach if i get it as there are multiple types of mutation. It would be more useful if i do pairwise comparison for specific cell types i guess

@StupidWolf for example I would like to see what set of genes which changes or are significant when i traverse from HSC to down the lineage tree such as CMP,GMP,Monocyte, but I have also included LSC(Leukemic stem cells) and Blast. But not sure if this is right comparison, i did normalize and took the data for clustering that is fine ,but hypothesis testing across heterogeneous group which are diseased with mutation would that be correct?

@StupidWolf this answer let me go through again https://bioinformatics.stackexchange.com/questions/13627/how-to-design-deseq2-lrt-model-with-individuals-nested-in-2-levels

Looking at your design with ATpoint's table, there's a major issue. some of your mutations have n=1 in some of the interactions, for example TET2-R550 or NRAS-G13D. or NPM1-L287ins(TCTG) is exclusive for blast. So it's hard for you to make these comparisons say HSC.. when it is so incomplete

I think you basically need to think about how to design this

"there's a major issue. some of your mutations have n=1 in some of the interactions, for example TET2-R550 or NRAS-G13D. or NPM1-L287ins(TCTG) is exclusive for blast. So it's hard for you to make these comparisons say HSC" now i get it at least i need a replicate ..may be i can make another categorical variable

> table(doldat$Group, coldat$Mutation)

        DNMT3A-R882H FLT3-ITD IDH2-R140Q KRAS-G13D MAU2-Q98H No NPM1-L287ins(TCTG) NRAS-G13D STAG2-R614 TET2-E1357 TET2-R550
  HSC              0        0          0         0         0  7                  0         0          0          0         0
  Blast            2        4          2         1         1  0                  4         1          1          6         1
  CMP              0        0          0         0         0  8                  0         0          0          0         0
  GMP              0        0          0         0         0  7                  0         0          0          0         0
  LSC              2        3          0         0         1  0                  0         1          1          0         0
  Mono             0        0          0         0         0  6                  0         0          0          0         0

Look at this table, you see by eye how some mutations are nested with Group. If you want to include it then you will need to subset the experiment into multiple separate runs, depending on the question you want to answer.

Lets say if i want to see the effect of mutation+Group  across HSC , LSC and Blast would the above metadata design work ?


text:
PCA analysis and nucleotide diversity pi
I have 25 individuals' whole-genome sequence (collected from 5 different locations) and after doing a variant calling analysis, I found 2734 SNP+INDELS marker (SNPs + INDELS = 2734 ). Then I plotted the PCA plot (pc1 and pc2) using variant calling data, i.e. raw SNPs.
Now I am trying to interpret the population structure.
Here I can see that these individuals from different locations are genetically different and there is no difference between individuals from the same geographic foci.
Are there any other things that I can report?
EDIT: Each circle(except green) represents a cluster of 4 individuals (that overlap with each other) and the green dot represent one individual

Then I calculate nucleotide diversity (π) values (across the whole genome) of each cluster observed in PCA plot: What is best way to show that information?
The mean pi value of each cluster:
1 5.61 E-05
2 4.76E-0.5
3 0.00045
4 2.56e-05
5 4.69e-0.5

Is there any test that should I perform to see if there are any significant differences, such as panmixia, migration or neutrality?
Can someone explain to me how to interpret and report the finding from PCA and π?

What's being shown on the PCA plot? I'd have expected to see the 25 samples, but it looks like you're plotting just the cluster centers, which makes it impossible to evaluate how tightly grouped the samples are. What's the meaning of the inner green circles?

HI, the sample that I plotted overlap each other that's why it is looking like that. Each circle represent (except green) a cluster of 4 individuals (that overlap) and the green dot represent one individual

@NuclearHoagie I made changes to the original post.

when you say you are interested in population structure, what do you mean? why are you estimating pi within clusters? what is the meaning of the clusters? did you place the clusters together or are they the output of some procedure? Is there a reason that you are not using $F_{ST}$ between clusters, if they are pre-defined? Have you considered using STRUCTURE or another tool that does unsupervised clustering (not just PCA projection)?

I'm happy to answer this if @Ashar is still seeking a solution


text:
What unit is the length of a protein measured in?
Sorry if this question is obvious - I'm not familiar with proteomics so I'm not sure of the standard units.
In this paper, it talks about using the length of a protein to calculate the NSAF value.  However, it doesn't specify the units.  Would it be peptides, amino acids or nucleic acids? Or something else?

Number of amino acids.

The paper clearly states protein length (number of amino acids).

Isn't this exactly the type of questions we like on this site?

As mentioned by others, it's the number of amino acids. Alternative measurements could be in Angstroms or similar, but that would require at least an approximate structure.


text:
Clusters in scRNA of Glioblastoma
There is a scRNA dataset on Glioblastoma from 10xgenomics.
Unfortunately, there seems to be no research paper related to that data.
Question 0 : May be there is paper  ? How to find it ? I am not familiar with 10xgenomics - may be there is some standard way to find papers related to their datasets ? May be their data also placed in GEO database, then what is GSE-id ?
Any way, we can see several clusters in the data (see figure above or switch to tab 'analysis' at the link above).
Question 1: What can be biological interpretation of these clusters ? I would expect cancer cells to be more homogeneous, but here we see several well-separated clusters, probably something like that has been already discussed in the literature ?
PS
That dataset has been placed on kaggle https://www.kaggle.com/reminho/human-glioblastoma-dataset
by enthusiast
and one can find several notebooks in R and Python to play with it. Clusters are also present there.
In general there are many datasets/notebooks/people related to bioinformatics on kaggle, it is quite convenient platform for everything related to data science - data sharing, data processing on their servers for free (like colab), commenting on the works of others and so on.
It seems it would be useful to create bioinformatics "community" on kaggle (there are already 12 other "communities") if you find that idea - useful please share your comments/make upvote/show yourself at:
https://www.kaggle.com/general/203136

I actually do not think cancers are homogeneous. There can be cancer cells, normal cells, various immune cells in cancer tissue. So it is not surprising to see these many clusters, it may be helpful to try to find markers for different clusters and identify the cell type.


text:
How does decrease the size of "grey" module in WGCNA framework?
Relate to finding soft power, I ran "pickSoftThreshold" function and the result of that(sft) is as follow:
      "Power"   "SFT.R.sq"  "slope"   "truncated.R.sq"  "mean.k."   "median.k."     "max.k."
"1"     1         0.30       0.92         0.89           2056.1      1903.5          4558.9
"2"     2         0.50      -1.52         0.91            548.1       445.5          1973.9
"3"     3         0.78      -2.05         0.96            188.3       128.3          1015.9
"4"     4         0.81      -2.21         0.97            77.21       42.3            594.9
"5"     5         0.82      -2.21         0.98             36.1       16.2            373.2
"6"     6         0.80      -2.18         0.96             18.8       6.9             245.8
"7"     7         0.78      -2.07         0.94             10.6       3.1             167.8
"8"     8         0.93      -1.75         0.99             6.4        1.5             121.7
"9"     9         0.94      -1.74         0.99             4.2        0.79            103.0
"10"    10        0.95      -1.69         0.99            2.87        0.42            88.2
"11"    12        0.97      -1.60         0.98            1.52        0.13            66.4
"12"    13        0.96      -1.56         0.98            1.17        0.07            58.7
"13"    14        0.96      -1.52         0.98            0.92        0.04            52.3
"14"    15        0.96      -1.48         0.98            0.75        0.02            46.9
"15"    16        0.95      -1.45         0.98            0.62        0.01            42.2
"16"    17        0.95      -1.41         0.98            0.53        0.01            38.1
"17"    18        0.95      -1.36         0.99            0.45        0.01            34.6
"18"    19        0.94      -1.33         0.99            0.40        0.00            31.5

My network type is 'Signed Hybrid' and cor type is 'bicor' as follows:
sft = pickSoftThreshold(datExpr, powerVector = powers, networkType = "signed hybrid", corFnc="bicor" ,corOptions = list(maxPOutliers =0.1), verbose = 5)

based on the sft result that I presented in this post and also mean connectivity plot, I selected power = 7 to construct a co-expression network for 5000 most connected genes and 58 colerectal cancer patients. Now, I have 5 modules as below:
green     brown    grey   turquoise    yellow
 550       345      3400       610        95

As you see, 3400 genes are in the "grey" module. My question is that are these genes just in the "grey" module is reasonable? if it's not reasonable, how can I modify my constructing process? do the selected power as 7 is correct?
I really appreciated it if anybody share his/her comment with me.

How homogene is your data? All of your 58 samples are colorectal cancer patients? Isn't there any subtype of colon cancer? Why do you filter your data to the 5000 more connected? The tutorials clearly say you shouldn't filter the data (Besides how did you select that?). Is your data normalized or preprocessed in any way?

Thanks@llrs. my 58 samples are as one subtype of the colorectal cancer patients. most connected genes are acceptable filter because Professor Horvath used this filter before. my data was normalized through TPM method. could you please share your comment with me?

That someone used it in some analysis doesn't mean it is right with your analysis. You didn't share how did you filter those genes (or [edit] them on the question), I do not have further advice without further information


text:
How to control/normalize for number of reads when calling SNPs using RNA-Seq?
I used the GATK pipeline to call SNPs on males and females using RNA-Seq data. But the males have a higher read count (~43-46M reads) than the females (~40-42M reads). This causes SNP counts to be higher for males and lower for females. This is strange because I am interested in X-chromosome and hence the question "How to control for the number of reads when calling SNPs using RNA-Seq?"

I think there is no way to do this, the ultimate solution is to get more data. If you have more reads, you naturally have more SNPs called. You can find the intersection of males and females if you want to analyze the two together.

would it be acceptable to subsample all the reads at 40M reads with the same seed and then do the GATK pipeline? I am thinking of seqtk to subsample.

That is not even 10% difference, plus you would need to check whether the one group has the same number of reads aligning to the exome. It could be that one group has more gDNA contamination, or rRNA, or whatever, and naive per-million scaling would fail. I am not sure though that all is necessary here, just run the pipeline RNA-seq is anyway suboptimal for variant calling and you will need to confirm interesting findings with an indepedent approach, so don't overthink it at this stage, run the pipeline and see whether you even find something interesting at all.

Some random thoughts:
A) How do you define a SNP in the X chromosome? I wouldn't expect to identify polymorphic positions in males, so I assume you call SNP a site in which the male individual is homozyogus for an allele different from the reference.
B) In line with what @ATpoint suggests, give a look at how many reads are aligned in males and females (especially in X chromosome).
C) Which organism is this? Are there any reasons to suspect that females are more different from the reference than males?


text:
Biology behind PCA analysis based on SNP
What is the biology behind SNP based PCA analysis to study population structure?  I am reading some articles where they used PCA analysis to compare isolates of drosophila that is collected from a different location. AND by PCA they found that SNP is different between all populations. What does that mean? (what it tell you about population structure/genetics?)

First of all, PCA is a technique for dimension reduction. Basically, the goal is to compare tens of thousands of SNPs in Drosophila. Now if you only have 2 SNPs, you can plot them on a 2D scatter plot. If you have 3 SNPs, you may try a 3D plot. But now imagine you have 30,000 SNPs, but you CANNOT plot a 30000-dimensional plot. To visualize this high dimensional data, what we can do is to perform dimensional reduction like PCA. PCA tries to find a set of orthogonal coordinations that explains most of the variation in the data (if there no variation, there is no information contained in the data, which essentially means there is no data). The idea is that PC1 carries most variation can be explained, and PC2 carries the second most... For lower PCs like PC50 or PC60, they probably only carry noise in the data. Therefore, the higher PCs (PC1, PC2 and so on) effectively summarizes the useful information in the data. So you can visualize the "structure" of the data in a 2D PCA plot.
By looking at the distance between points on a PCA plot, you can tell how similar the two data points are. But if you see two populations that are perfectly separated on PCA plot, it does not mean that the 2 population differ completely at every SNP, because PCA is a summarization of all SNP included.

I understood that PCA is performed to reduce dimensions. I didn't get  ""Now if you only have 2 SNP....."". let say I got 2 SNP and I want to plot that on scatter plot. What will the x and y-axis be?  And why with every SNPS dimensions are increasing. As you said that for 3 SNPs, I can use a 3d plot. What is the relation between SNPs and dimension?

one more question: Let say there are two population of drosophila A and B collected from 2 different locations. Both populations contain 5 isolates each.  After performing Variant calling analysis (for both populations at once) I got a VCF file showing that I have 5000snps. Here how would I know how many SNPs belong to each population? ( I am assuming 5000snps is detected from both 2 population collectively)

@Ashar please try not to ask additional questions in the comments. Add on any minor points to the main question or create a new question. This makes it easier for other people to search for questions. Thanks.

SNPs don't belong to one population. Most of them differ by allele frequency between two population, unless one SNP is fixed in one population (which is very rare).


text:
ggplot2 geom_bar() , do not merge y-values shared same x-ids
I am plotting a barplot using ggplot + geom_bar() function in R.
I have some the recoreds which share the same x-axis identifiers.
geom_bar() did merge all values that shared the same x-axis identifiers.
Like the following example.
library("ggplot2")

df <- data.frame(x = c("a","a","b","c","d","e"),
                 y = c(5,5,5,5,5,5),
                 stringsAsFactors = FALSE)

ggplot(df, aes(x, y)) + geom_bar(stat = "identity")

The erroneous graph is displayed below.
The highest bar means the sum of multiple y-values shared same x-axis names.
Anyone else could help me to solve this problem? Thanks.

Please provide a reproducible example so we can understand what's going on.

Interesting question. Your data is out of shape though - you cannot plot duplicate data points on an axis - how will anyone reading the plot know the difference between the two data points?

Finally， I solved it using barplot() function that does not merge values that have same x-axis name.

As you seem to have surmised, the graph is not erroneous, but correct. You have a total of 10 for a and that's what it's showing. If you want it to ignore duplicate X-axis values then you need to tell ggplot which datapoints it should actually be using. barplot() is assuming that your duplicate a values are intentional, ggplot assumes they should be grouped. These are both reasonable assumptions and it's your responsibility to determine beforehand which is correct.


text:
OpenAPI-Generator Python Client Instructions Examples
Trying to learn a little bit of Python and Biopython I ended up figuring out what a Rest API is when the RCSB PDB repository released its new Rest API.
I tried to retrieve info from the PDB repository using python3 requests,
then I managed to use OpenAPI-generator to generate what they call a 'Python3 Client'.
Now I am at loss I have this ApiClient class:
print dir():

 ['NATIVE_TYPES_MAPPING', 'PRIMITIVE_TYPES', '_ApiClient__call_api',
  '_ApiClient__deserialize', '_ApiClient__deserialize_date',
  '_ApiClient__deserialize_datetime', '_ApiClient__deserialize_file',
  '_ApiClient__deserialize_model', '_ApiClient__deserialize_object',
  '_ApiClient__deserialize_primitive', '__class__', '__delattr__',
  '__dict__', '__dir__', '__doc__', '__enter__', '__eq__', '__exit__',
  '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__',
  '__init__', '__init_subclass__', '__le__', '__lt__', '__module__',
  '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__',
  '__setattr__', '__sizeof__', '__str__', '__subclasshook__',
  '__weakref__', '_pool', 'call_api', 'client_side_validation',
  'close', 'configuration', 'cookie', 'default_headers',
  'deserialize', 'files_parameters', 'parameters_to_tuples',
  'pool', 'pool_threads', 'request', 'rest_client',
  'sanitize_for_serialization', 'select_header_accept',
  'select_header_content_type', 'set_default_header',
  'update_params_for_auth', 'user_agent']

and don't know how to use it. Googling for examples or searching for info on the OpenAPI-generator website doesn't seem to work for me, I find plenty of tutorials and examples on how to generate the client but none about how to use it.
The generated client class should help to better interact with the server as I understood (I could be wrong here). I am asking for helps and hints here because the OpenAPI webpage under Community points to stackoverflow under [openapi-generator] tag. Any tutorial, easy documentation that could help me understand the output of the generator?

OK fine, found useful information on READme.md generated by OpenAPI-generator itself
and tried pdoc to get html document of my OpenAPI-generator generated client Module
Moreover on stackoverflow How to generate API Documentation using openapi-generator? you can find how OpenApi Generator can do generation of html documentation from your yaml or json definition of your API.
For openAPI json api specification I found RapiPDF too : Web Component to generate PDF from OpenAPI Spec (json yml file).


text:
How to filter a genome assembly consistsing of a large number of contigs?
I did some de novo genome assemblies with Illumina PE data using SPAdes, whereas most of them consisting of a large number of contigs(>1000). I have several questions below.

Do we need to filter short contigs?
How to determine the appropriate length threshold？

Hi, YP CHEN, could I know what kind of genome it is? bacterial? viral? eukaryotic? How long is the genome?

Fungal genomes， genome sizes are about ~ 50M

Thanks. I can't help in that case. But I guess what determines your length threshold is what use you intend to make from that assembly.


text:
How to select specific lines and specific columns
I have  a data frame that contain lot of columns and lines, like this for example
Column1   Column2
5            T1
1            T2
1            T3
1            T1
2            T1
1            T1
2            T1
1            T1

I want to extract just lines that have the value 1 in column 1 and Value T1 in Column 2
So the result I want is
Column1   Column2
1            T1
1            T1
1            T1

Can you tell me please how I can do that

Hi @James what code have you written so far? You must have imported into R and done some basic dataframe work.

Hello , thank you very much for you answer , yes i tried this one : Pokemon_2 <- dropRows(ds, Pokemon$Generation == 1 , Pokemon$Type_1=="Water" )   but it didn't work  , Error in dropRows(ds, Pokemon$Generation == 1, Pokemon$Type_1 == "Water") :
could not find function "dropRows"

You should look at the subset function. Please google "filter data frame by column value" - much like writing for loops in basic programming languages, it is a basic operation that does not warrant creating new posts on forums.

I’m voting to close this question because it is a basic R operation has nothing to do with bioinformatics specifically

You can simply use subset function in R and use the code given below,
    new_dataframe <-
    subset(name_of_your_dataframe,Column1 == 1 &&
    Column2 ==  "T1")

Please keep in mind that R is case sensitive and you should be careful while using the variable names.


text:
What is the best way to start for analysing DNA methylation TGCA datasets using R? I want to learn but I'm a novice
I am relatively new to programming and only have a basic knowledge base for R. I am keen on learning how to manipulate DNA methylation arrays for analysis related to Glioblastoma Multiforme. More specficially I will want to try and find patterns and trends associated with differences in survival linked to expression of mutations and how significant the relationship is with progression free and overall survival.
Any help would be MASSIVELY appreciated. Thank you!

Does the R library you want to use have a vignette?