text:
ClustalOmega alignment
I am performing  multiple sequence alignment using Clustal Omega for 600 RNA covid-sequences, i.e. cDNA on Genbank, with input characters around 30574 characters for each sequence.
I am running it on windows cmd.
I defined parameters in input of maxseqle = 37000, however, the output each time gives duplicated length around 75000 charcters per each sequence.
clustalo.exe -i Allseq.fasta --is-profile --use-kimura --seqtype DNA --maxseqlen 37000 --threads 8 -o myclustalv3.fasta

How can I solve this problem, what will I define in input parameters. This problem occurs with large no of sequences 600 sequences. In test of 40 seq, it gave me normal output length around 32000 characters but when align all 600 sequences together, it gave this duplicated length results.

Same question here https://www.researchgate.net/post/Problem_in_Clustal-Omega_alignment

HEre's your response on Researchgate: "I concern the length because I am doing comparative analysis betweem Muscle,Mafft and ClustalOmega and the length of clustal is totally different than the others which made the analysis of alignment is totally different"

Have you had a look at your alignment? I personally use Jalview but any alignment viewer should work. Colour > Percentage Identity (or Nucleotide) should help with visualization.
Does the alignment look reasonable (and gaps are just inserted at the end for some reason), or is it pretty messy? If it looks messy, is it just a few outlier sequences that cause the problem, or does it look like you have two main groups of sequences that align between themselves but not with the other? In the latter case, could it be possible that some of your initial sequences were inputed in the reverse orientation? I'm not sure about the command line implementation of Clustal Omega, but I know that at least the versions I use don't automatically reverse-complement sequences to be aligned.

This is the right answer, I personally would avoid Jalview but its a matter of opinion.

Firstly its good you've done some proofing of your output wherein you spotted,

75000

The alignment output is approximately double the genome size, just over double the input.
An off-question point, I assume the computer you are working on has at least 8-cores because thats what you are requesting.
Clustal omega is a very solid algorhithm and has been around for a long time in its current omega version. The error you are encountering is an input sequence bug. @Laura's suggest is correct and is formal practice: you must manually verify bioinformatics output and in this case an alignment viewer is the correct proceedure.
Thereafter I suspect you would find the answer, such as omitting carriage returns in the fasta format, seeing taxa names caught up post 30K.
If you run,
wc ./myclustalv3.fasta # output file
wc ./Allseq.fasta  # input file
wc ./testgenome.fasta

You will probably find that the number lines per file does not agree between the input and output format OR the number of lines of both input and output is <<600 or <<1200 (600 x2 ), i.e. the expected number of sequences. This points to parsing errors within the input data. Fasta has a line for the sequence identifier and a line for the sequence. It is using > for the sequence identifier, but will use a carriage return for sequence.
Again manually assessing an alignment is baseline proceedure regardless of output.
An example (guess) about what might have happened is that every other input has omitted ">", omega will then see the first ">" then read everything until it sees another ">" as sequence data. Thus its doubled the sequence length and dragged a sequence identifier exactly in the centre of the doubled up genomes.

An issue with fasta is there is no internal proofing of the format. In other formats the number of sequences and their length is specified within the file, so it is easy to internally verfiy the data prior a calculation. Code-wise it is known as 'setting' and you can't do that with fasta.


text:
Help to create a dataframe in Python from a FASTA file
I want to create a dataframe in Python starting from a FASTA format file. Given the toy FASTA file that I am attaching, I built this program in Python that returns four colums corresponding to id, sequence length, sequence, animal name and rows corresponding to all the data available.
However, I am trying to understand how to modify this code in order to create a dataframe in which classes Human and Dog have the same number of data. For example, I want to say to Python: "Append to record (that is the empty list) id, sequence length, sequence and animal for Human, but do it a number of times that is equal to the number of data of the class with minimum number of data (that is Dog)".
I think that a while loop is needed but I have a bit troubles to understand how to do it. Any suggestion ?
Below the Python code I wrote and the FASTA format file I used.
import pandas as pd
import re
def read_fasta(file_path, columns) :
    from Bio.SeqIO.FastaIO import SimpleFastaParser
    with open("Proof.txt") as fasta_file :
        records = [] # create empty list
        for title, sequence in SimpleFastaParser(fasta_file): #SimpleFastaParser Iterate over Fasta records as string tuples. For each record a tuple of two strings is returned, the FASTA title line (without the leading ‘>’ character), and the sequence (with any whitespace removed).
            record = []
            title_splits=re.findall(r"[\w']+", title) # Data cleaning is needed

            record.append(title_splits[0])  #First values are ID (Append adds element to a list)
            record.append(len(sequence)) #Second values are sequences lengths
            sequence = " ".join(sequence) #It converts into one line
            record.append(sequence)#Third values are sequences

            #Fourth column will contain the species
            if "Human" in title_splits:
                    record.append("Human")
            else:
                    record.append("Dog")

            records.append(record)
    return pd.DataFrame(records, columns = columns) #We have created a function that returns a dataframe

#Now let's use this function by inserting in the first argument the file name (or file path if your working directory is different from where the fasta file is)
#And in the second one the names of columns
data = read_fasta("Proof.txt", columns=["id","sequence_length", "sequence", "animal"])
data

The FASTA format file is this:
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Dog|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Dog|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Dog|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Dog|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK
>Numer|Human|Hearth
HSSFIEIVNIEHVIEHIVK

My code prints a dataframe like:
       id  sequence_length                               sequence animal
0   Numer               19  H S S F I E I V N I E H V I E H I V K  Human
1   Numer               19  H S S F I E I V N I E H V I E H I V K  Human
2   Numer               19  H S S F I E I V N I E H V I E H I V K  Human
3   Numer               19  H S S F I E I V N I E H V I E H I V K  Human
4   Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
5   Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
6   Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
7   Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
8   Numer               19  H S S F I E I V N I E H V I E H I V K  Human
9   Numer               19  H S S F I E I V N I E H V I E H I V K  Human
10  Numer               19  H S S F I E I V N I E H V I E H I V K  Human
11  Numer               19  H S S F I E I V N I E H V I E H I V K  Human

But I would like that the number of rows for Human is the same for Dog (because, in other words, I would like the same number of data for each class).
Hoping to have been clear, I thank you in advance.
P.S.: My program is an adapted version of this notebook .

Please correct me if I am wrong, you want to print the lines Human 4 times because it is equal to the number of rows that corresponds to Dog, is that right?

Assuming what I said in the comments is true. Here's my attempt, there may be an easier solution but this should do it:
import pandas as pd
import re

# a function to return the count number of animals in the fasta in a dictionary
def count_animals(file_path, separator, animal_index):
  animal_dictionary={}
  with open("Proof.txt") as fasta_file:
    fasta_contents = fasta_file.readlines()
    for line in fasta_contents:
      if line.startswith(">"):
        animal = line.split(separator)[animal_index]
        if animal not in animal_dictionary:
          animal_dictionary[animal] = 1
        else:
          animal_dictionary[animal] = animal_dictionary[animal] + 1
  return animal_dictionary

def read_fasta(file_path, columns) :
  counter = 0  # counter used to print no more than number of times 'Dog'
  animal_d = count_animals(file_path, '|', 1) # execute function
  max_animal =  max(animal_d.values()) # get max number of animals
  min_animal =  min(animal_d.values()) # get min number of animals

  from Bio.SeqIO.FastaIO import SimpleFastaParser
  with open("Proof.txt") as fasta_file :
    records = [] # create empty list
    for title, sequence in SimpleFastaParser(fasta_file): #SimpleFastaParser Iterate over Fasta records as string tuples. For each record a tuple of two strings is returned, the FASTA title line (without the leading ‘>’ character), and the sequence (with any whitespace removed).
      record = []
      title_splits=re.findall(r"[\w']+", title) # Data cleaning is needed

      record.append(title_splits[0])  #First values are ID (Append adds element to a list)
      record.append(len(sequence)) #Second values are sequences lengths
      sequence = " ".join(sequence) #It converts into one line
      record.append(sequence)#Third values are sequences
      # if 'Human' exists in title and the counter is less than min_animal('Dog'), add the record
      if 'Human' in title and counter < min_animal:
        counter += 1 # add +1 to counter

        #Fourth column will contain the species
        record.append("Human")

        records.append(record)
      elif 'Dog' in title:
        record.append("Dog")

        records.append(record)

    # print(records)
  return pd.DataFrame(records, columns = columns) #We have created a function that returns a dataframe

#Now let's use this function by inserting in the first argument the file name (or file path if your working directory is different from where the fasta file is)
#And in the second one the names of columns
data = read_fasta("Proof.txt", columns=["id","sequence_length", "sequence", "animal"])

print(data)

Should give:
      id  sequence_length                               sequence animal
0  Numer               19  H S S F I E I V N I E H V I E H I V K  Human
1  Numer               19  H S S F I E I V N I E H V I E H I V K  Human
2  Numer               19  H S S F I E I V N I E H V I E H I V K  Human
3  Numer               19  H S S F I E I V N I E H V I E H I V K  Human
4  Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
5  Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
6  Numer               19  H S S F I E I V N I E H V I E H I V K    Dog
7  Numer               19  H S S F I E I V N I E H V I E H I V K    Dog

Thank you very much !

An example of contructing a pandas dataframe for fasta on Stack Exchange is here . This is clean code deploying zip outside Biopython. Constructing a pandas dataframe directly via dictionary.
Are you intending to use split method directly on the sequence, such that each residue is a separate column in a dataframe?
Anyway,
title_splits=re.findall(r"[\w']+", title) # Data cleaning is needed

This is a bad idea because you will remove the Genbank unique identifier (its an integer) and this should not leave the sequence.
You will need to switch out of the Fasta object. Your approacha and @user324810's use of
sequence = " ".join(sequence)

is okayish. However, I think it should be '', or "" (no space character). I spotted manually removing fasta, well the method should do this automatically. I am not familiar with Fastaparser, but AlignIO and SeqIO will both automatically process fasta without the need deal with ">". Biopython will have a length method and whilst you have the Biopython object it is preferable to use Biopython (OOP) methods.
from Bio.SeqIO.FastaIO import SimpleFastaParser

This should be at the top of the code.
Overall, you manually open the data, pass it through the Bioparser then immediately dump the object. I do agree to write a pandas dataframe from a Biopython object these conversion to string, but Biopython will provide a method to print the object (pass it to a string) without all the parsing. Some of the parsing is removing key data (Genbank ID).
Personally I would look at the Stackoverflow example and base your code around this.
You clearly have an understanding of pandas, but I get the impression Biopython is something your not familiar with. My suspicion is that Biopython would infact answer the entirety of your question (whatever it might be) without passing to pandas.
pandas is cool and powerful for non-molecular data, but unless your using machine/deep learning (where a dataframe is necessary input) I think Biopython will provide most of the answers.
Important Terminology
The class is FastaIO and the method is SimpleFastaParser. Humans and dogs are instances of this method and are therefore object(s).
Please don't write 'class' of dogs and humans in a formal report because it is incorrect OOP and will cause alarm if your reviewer is a coder. Your code is functional programming rather than OOP. You could say the human/dog fasta object, or just refer to it as human/dog dataframes and avoid any OOP terminology.

Thank you for the suggestions @M__ but SimpleFastaParser does not remove Genbank ID... or at least the FASTA file that I have is with the header like this: >Numer|Animal|OtherInformations . And once SimpleFastaParser, for each record, has returned a tuple of two strings (and the header that is without ">"), re.findall has splitted the header string according to all characters (or at least according to the "|") that I am interested in, I can use this splitting to identify the name of the species and append it to the empty list...

The only thing is in fact that re.findall splits also according to point as separator and I noticed that GenBank IDs are numbers like NP288364.1 ... So re.findall will split NP288364.1 into NP288364 and 1 . Is it so terrible ? I found that the ".1" indicates the "version" of the sequence. Do you know a way to change re.findall so that I can split without considering the points ?

Yes Biopython is new for me but it does not delete the ID... maybe the fault lies with re.findall, isn't it ?

Yes, you should use different parsing such as re.sub explicitly stating "|" such as using a capture "(.+)|.*" and $1, or \1 in re.sub or using or re.final ([\w|\d], string).

Just to add, usually the way to move a sequence object to a neat, decluttered string in Biopython is not to do this manually as done here:
title_splits=re.findall(r"[\w']+", title)
sequence = " ".join(sequence)

both being from SimpleFastaParser(fastafile)
It should be possible and, if so, better to state,
 str(title)
 str(sequence)

If you can issue a print command and the string appears nice and neat on the screen str() will work.
This becomes important for complex sequence/alignment objects, such as SeqIO and AlignIO. The reason it works is __str__ has been declared or available through OOP inheritance.

Thank you, however in fact I realize only now that sequence = " ".join(sequence) is not necessary since SimpleFastaParser already joins sequences into one line and as regards the header I tried with str as you say but it does not work because it leaves it as it is Numer|Human|Hearth and I do not know how to say "append Number". I could write like  title_str=str(title) record.append(title_str[0:6]) but the length of the id number varies for each header...

Fine good to know. str isn't implemented for this given method (they haven't written one). It is a standard method in most OOP, but not all

source: Stack Exchange

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

sort: votes


text:
How to process multiple VCF files present in a 2 folders and create afd plots
This question has also been asked on Biostars
I have 2 directories that contains multiple VCF files.

I want to write code in python that will process all vcf files and create an AFD (allele frequency difference) plot for each VCF file present in a folder
I have started with this code, which attempts to convert the VCF files into a pandas dataframe:
# Import Module
import os
import pandas as pd
# Folder Path
path1 = "C://Users//USER//Desktop//isiah/VCFs_1/"
path2 = "C://Users//USER//Desktop//isiah/VCFs_2/"
#os.chdir(path1)

# Read text File

def read_text_file(file_path1,file_path2):
    with open(file_path, 'r') as f:
        print(f.read())

# iterate through all file
for file in os.listdir():
    # Check whether file is in text format or not
    if file.endswith(".vcf"):
        file_path1 = f"{path1}\{file}"
        file_path2 = f"{path2}\{file}"
        print(file_path1,"\n\n",file_path2)

# call read text file function
#data = read_text_file(path1,path2)
#print(read_text_file(path1,path2))
Final= pd.DataFrame(read_text_file(path1,path2))
print(Final)

and its giving me output without calling function and also gives me empty dataframe. Also i am not sure if result is showing all vcf files record or not.
Kindly help me to modify this code. I need my all result in dataframe so i can make plot easily
Folder 1 contains 3 test vcfs whereas folder 2 contains 4 test vcfs:

This is test 1 vcf file:
##fileformat=VCFv4.2
##ALT=<ID=NON_REF,Description="Represents any possible alternative allele at this location">
##FILTER=<ID=FS,Description="FS>60.0">
##FILTER=<ID=HS,Description="HaplotypeScore<13.0">
##FILTER=<ID=LowQual,Description="Low quality">
##FILTER=<ID=MQ,Description="MQ<40.0">
##FILTER=<ID=MQRS,Description="MQRankSum < -12.5">
##FILTER=<ID=Mask,Description="Doesn't overlap a user-input mask">
##FILTER=<ID=PASS,Description="All filters passed">
##FILTER=<ID=QD,Description="QD<2.0">
##FILTER=<ID=RPRS,Description="ReadPosRankSum < -8.0">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=PGT,Number=1,Type=String,Description="Physical phasing haplotype information, describing how the alternate alleles are phased in relation to one another">
##FORMAT=<ID=PID,Number=1,Type=String,Description="Physical phasing ID information, where each unique ID within a given sample (but not across samples) connects records within a phasing group">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
##FORMAT=<ID=RGQ,Number=1,Type=Integer,Description="Unconditional reference genotype confidence, encoded as a phred quality -10*log10 p(genotype call is wrong)">
##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias.">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##INFO=<ID=BaseQRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt Vs. Ref base qualities">
##INFO=<ID=ClippingRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref number of hard clipped bases">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth; some reads may have been filtered">
##INFO=<ID=DS,Number=0,Type=Flag,Description="Were any of the samples downsampled?">
##INFO=<ID=ExcessHet,Number=1,Type=Float,Description="Phred-scaled p-value for exact test of excess heterozygosity">
##INFO=<ID=FS,Number=1,Type=Float,Description="Phred-scaled p-value using Fisher's exact test to detect strand bias">
##INFO=<ID=HaplotypeScore,Number=1,Type=Float,Description="Consistency of the site with at most two segregating haplotypes">
##INFO=<ID=InbreedingCoeff,Number=1,Type=Float,Description="Inbreeding coefficient as estimated from the genotype likelihoods per-sample when compared against the Hardy-Weinberg expectation">
##INFO=<ID=MLEAC,Number=A,Type=Integer,Description="Maximum likelihood expectation (MLE) for the allele counts (not necessarily the same as the AC), for each ALT allele, in the same order as listed">
##INFO=<ID=MLEAF,Number=A,Type=Float,Description="Maximum likelihood expectation (MLE) for the allele frequency (not necessarily the same as the AF), for each ALT allele, in the same order as listed">
##INFO=<ID=MQ,Number=1,Type=Float,Description="RMS Mapping Quality">
##INFO=<ID=MQRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities">
##INFO=<ID=QD,Number=1,Type=Float,Description="Variant Confidence/Quality by Depth">
##INFO=<ID=RAW_MQ,Number=1,Type=Float,Description="Raw data for RMS Mapping Quality">
##INFO=<ID=ReadPosRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias">
##INFO=<ID=SOR,Number=1,Type=Float,Description="Symmetric Odds Ratio of 2x2 contingency table to detect strand bias">
##reference=file:///Users/sian_bray/Dropbox/Yant/Cochlearia/02_Data/References/pseudohap_cochlearia_135M_3kb.fasta
##source=SelectVariants
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  ALO_006 ALO_007 ALO_013 ALO_017 BEA_002 BEA_004 BEA_010 BRI_002 BRI_005 BRI_006 BRI_009 ELI_001 ELI_002 ELI_003 ELI_004 FRE_013 JON_001 KVA_002 KVA_003 KVA_009 KVA_010 LAB_004 LAB_300 LAB_400 LAB_500 LNL_001 LNL_002 LNL_003 LNL_008 MEL_001 MEL_002 NEN_001 NEN_003 NEN_200 NEN_300 ROT_004 ROT_006 ROT_007 ROT_013 RUZ_001 SAL_003 SAL_004 SAL_006 SKF_002 SKF_003 SKF_005 SKF_009 SKI_004 SKI_005 SKI_008 SKN_001 SKN_002 SKN_005 SKN_008 SPU_006 SPU_008 SPU_009 SPU_010 TET_002 TET_004 TET_006 TET_008 TRO_001 TRO_003 TRO_005 TRO_009 VAG_003 VAG_004 VAG_007 VAG_102 VEG_003 VEG_004 WOL_002 WOL_006 WOL_009 WOL_010
237 50041   .   G   T   4099.26 PASS    AC=35;AF=0.133;AN=264;BaseQRankSum=0.00;ClippingRankSum=0.00;DP=1228;ExcessHet=4.3562;FS=0.644;InbreedingCoeff=-0.0953;MQ=57.84;MQRankSum=0.00;QD=25.94;ReadPosRankSum=-1.570e-01;SOR=0.768 GT:AD:DP:GQ:PL  0/0/0/0:11,0:11:14:0,14,33,66,425   0/0/0/0:11,1:12:0:0,0,0,27,369  0/0/0/0:13,3:16:0:0,0,0,0,355   0/0/0/0:41,0:41:50:0,50,120,241,1800    0/0/1/1:5,8:13:1:293,15,0,1,170 0/0/0/1:5,3:8:0:104,0,0,10,191  0/0/0/0:8,0:8:10:0,10,24,48,272 0/1/1/1:1,8:9:11:318,33,11,0,26 ./././.:0,0 1/1/1/1:0,2:2:2:87,12,6,2,0 0/0/0/0:5,2:7:0:0,0,0,0,120 0/0/0/0:16,2:18:0:0,0,0,21,515  0/0/0/0:8,1:9:0:0,0,0,13,270    0/0/0/0:15,3:18:0:0,0,0,0,459   0/0/0/0:10,1:11:0:0,0,0,21,363  0/0/0/0/0/0:15,0:15:12:0,12,26,45,71,116,589    0/0/0/0/0/0:7,2:9:0:0,0,0,0,0,0,151 ./././.:0,0 0/0/0/0:14,0:14:14:0,14,33,66,495   0/0/0/0:1,0:1:1:0,1,3,6,42  0/0/0/0:9,0:9:10:0,10,24,48,360 0/0:21,2:23:14:0,14,726 0/0:14,2:16:11:0,11,452 0/0:14,2:16:0:0,0,467   0/0:23,9:32:0:0,0,527   0/0/0/0:10,0:10:12:0,12,30,60,359   0/0/0/0:11,1:12:0:0,0,0,27,404  0/0/0/0:13,0:13:16:0,16,39,78,526   0/0/0/0:20,3:23:0:0,0,0,2,655   1/1/1/1:0,4:4:5:155,24,12,5,0   0/0/0/0:2,3:5:0:0,0,0,0,8   0/0:23,2:25:0:0,0,798   0/0:12,2:14:0:0,0,381   0/0:14,0:14:36:0,36,540 0/0:24,0:24:63:0,63,945 0/0/0/0:8,1:9:0:0,0,0,12,245    0/0/0/0:27,0:27:32:0,32,78,157,1170 0/0/0/0:18,5:23:0:0,0,0,0,474   0/0/0/0:22,1:23:15:0,15,50,114,736  0/0:10,0:10:27:0,27,405 0/0:8,0:8:24:0,24,322   0/0:13,0:13:36:0,36,540 0/0:4,0:4:12:0,12,155   0/0/0/0/0/0:5,2:7:0:0,0,0,0,0,0,125 0/0/0/0/0/0:19,1:20:0:0,0,0,18,50,106,731   0/0/0/0/0/0:28,3:31:0:0,0,0,0,1,82,978  0/0/0/0/0/0:53,10:63:0:0,0,0,0,0,0,1843 0/0/0/0:1,1:2:0:0,0,0,0,4   1/1/1/1:0,22:22:27:917,132,66,27,0  ./././.:0,0 0/0/0/0:15,0:15:17:0,17,42,84,630   0/0/0/0:14,2:16:0:0,0,0,27,421  0/0/0/0:17,0:17:21:0,21,51,102,670  0/0/0/0:17,7:24:0:0,0,0,0,323   0/0:11,2:13:0:0,0,380

Answer from @user324810, converted from comment:
There is an alternative than doing it with Python unless it is absolutely necessary. If you want to merge multiple samples into one, use bcftools merge. If you want simply to concatenate multiple VCFs, use bcftools concatenate. Once you obtain the output file, you can load the dataframe and plot it as you wish.

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


text:
Does it make sense to compare the mean relative coverages for contigs / MAGs from samples that were not co-assembled?
I have a subset of high-quality MAGs from different environments, sequenced on different runs. I want to know, for any particular MAG, if it (or something similar to it) is present in any particular sample.
I've recently started learning anvi'o for metagenomic visualization. One of its perks is that it can nicely represent coverages for either contigs or bins (tutorial example: https://merenlab.org/images/miscellaneous/2017-05-11-anvi-refine-by-veronika/Image_1.png , the circles marked Sample 07 and 08).

The coverage options are defined as follows:
Mean coverage
Average depth of coverage across contig.
Add up the coverage of each nucleotide in a contig, and divide by the length of the contig. This is the default view in most cases.

Mean coverage Q2Q3
Average depth of coverage across contig excluding nucleotide positions with coverage values falling outside of the interquartile range for that contig.
Calculated the same as mean coverage, except only incorporating those nucleotide coverage values that fall within 2nd and 3rd quartiles of the distribution of nucleotide coverages for that contig. This can help smooth out the mean coverage visualization by removing nucleotide coverage values from the equation that may be outliers due to non-specific mapping.

Additional representation options:
Abundance
Mean coverage of a contig divided by that sample’s overall mean coverage.
Abundance for a contig is constrained to within one sample. In a sense this view is telling you that those contigs with larger abundance values are more represented in that sample (i.e. recruited more reads) than those contigs with smaller abundance values. And it does this by providing the ratio of that contig’s mean coverage to that sample’s overall mean coverage incorporating all contigs. So if your abundance value is 2, that contig’s mean coverage is twice that of the mean for all contigs in that sample.

Relative abundance
Proportion of reads recruited to a contig out of the total reads recruited to that contig across all samples.
Relative abundance considers one contig across all samples. This will tell you in which sample a particular contig recruited the most reads. Since it is normalized to the total number of reads recruited to contig across all samples, values from all samples for a given contig will always sum to 1.

Max-normalized ratio
Proportion of reads recruited to a contig out of the maximum number of reads recruited to that contig in any sample.

This also considers one contig across all samples. But in this case the value is normalized to the single maximum value for that contig (as opposed to the sum as in relative abundance). Because of this the sample containing the contig that contributed the max value will always equal 1, and the value for that contig in the other samples will be the fraction of that max.
I'm trying to cross-map the reads between samples to the best MAGs (+within the same sample). The reason I'm trying to map to MAGs as opposed to genomes, as usual, is that some of my MAGs of interest represent genera (at the least) without available genomes, so it feels like I'd lose considerable information that way. But then, I'm learning this as I go along, so for all I know it might be a bad idea.
From what I understood, these settings generally make sense / are used for samples that have been co-assembled. However, if I've been assembling samples independently and co-binning (due to wanting to avoid chimeras) or I want to compare some of  the metagenomic composition of different environments (which really wouldn't make sense to either co-assemble or co-bin), would it still make sense to use coverages as a proxy for presence / absence / abundance?
Thank you for your time!

any normalization? 2. have you cross-mapped reads between samples? Basically, for metagenomics, if you want to compare abundance between samples, you usually have a huge genome database that you map different samples' reads to, and then you use that to tell the difference between samples in a more or less unbiased way. Is that what you are using anvi'o to do?

@MaximilianPress 1) I don't know. Would something like 'rel. Abundance' (added above) count? 4) Yes/no.


text:
How to get a list of all human morphogens?
Gene Ontology provides many genes annotated as taking part in various morphogenesis processes but I want to get a list of all morphogen coding genes specifically. Uniprot  does not have "morphogen" keyword.
I need some gene annotation database that has genes labeled as "morphogen".


text:
PDBQT file format
I am trying to understand the source code of the AutoDock Vina, but as I am new to the field of Bioinformatics it is a little hard to understand. It seems like they are using a tree data structure to represent the ligand and the receptor. I cannot understand what nodes and branches are in this data structure (Heterotree). I think in order to understand it I should first understand what information the PDBQT file contains. Can you tell me what branch and root is in PDBQT file format?

For these very specific source code questions, it makes a lot easier if you link the code you talk about as specifically as you can (i.e. as I did in the edit. GitHub allows you to link individual files). Anyway, good luck :-)

SE is a not a forum for personal opinions, but I feel like I need break the rules here... We see occasionally CS folk ask questions here wherein they are in a wet lab biochem group and are given a computational biochemistry task with little technical supervision and they frequently flounder over and over at the basics. if this premise is correct, I would strongly suggest seeking a relevant collaborator within your institute/department.

what information the PDBQT file contains

The PDBQT format is described here: https://autodock.scripps.edu/wp-content/uploads/sites/56/2021/10/AutoDock4.2.6_UserGuide.pdf —ignore the expired SSL certificate warning as I am pretty sure the Scripps Institute is not a criminal organisation.

what branch and root is in PDBQT file format

It is rather self evident but with a twist and as you have guessed it has to do with trees in graph networks. Lines of text <-> graph networks describing a molecule is a common nuisance in compchemistry. Normally the aim is connectivity, but here it is torisons. But the former is a nice example, a small molecules can be represented with SMILES (strings), where say isopropanol is CC(O)C, in which you have a branch going off in round brackets, whereas everything else is assumed connected sequentially. In the case of the PDBQT format you are describing what atom subsets of the molecule rotate together and actually not how they connect as given away by the fact that you don't have a special syntax for ring closures (in a SMILES you write cyclopropane as C1CC1, where the third carbon connects to C1, which was declared first. This is because the topology is dealt with elsewhere.
Resources
AutoDock 4 and AutoDock Vina are divergent pieces of software, but they have several common elements like the PDBQT format, so you get double the documentation (and papers to read)!


text:
What does it mean by "PC" in this article?
What does it mean by "PC" in this article?

a fast and efficient Cα-based method for accurately assigning protein secondary structure elements

Answer from @user438383, @matteo-ferla and @m, converted from comments [anyone, feel free to edit this answer to improve it]:
According to the intro it means 'pseudocenter'. All measurements are relative to an assigned point because a protein doesn't have a 'centre' they assigned it as a pseudocentre.
Note that PC normally means principal component, but that does not seem to be the case here.


text:
What is the cutoff distance in this algorithm?
In the following article:

a fast and efficient Cα-based method for accurately assigning protein secondary structure elements

What are the minimum and maximum cutoff distances in this algorithm?
Also, what does it mean by:

I.e. where does (2x43)x3 come from?

The 43 features are listed in the supplementary materials on page 4. For example,
Feature Type Atoms Involved
1 Distance i, i-5
2 Distance i, i-4
3 Distance i, i-3
4 Distance i, i-2
...
11 Angle i-1, i, i+1
12 Dihedral i, i+1, i+2, i+3

Are some of them; so Feature 11 is the angle between atom(i-1), atom(i), and atom(i+1). Since it says these were calculated for both the C-alphas and pseudo-centers, this gives 2 * 43. Then this is multiplied by three as mentioned in your quoted paragraph.
These features (measurements, almost) are numbers, not booleans - so I do not think that they use distance cutoffs as such, although the random forest method is presumably learning various cutoffs and encodes those in the structure of the forest.


text:
Find the genomic locations of exons of certain human gene
So I know Consensus CDS has the genomic locations of exons, however, only of the exons in the CDS. Everything before or after is not included.
NIH variation viewer has those (kinda-) radio buttons next to the gene name, and when you hover over them, you do get to see the genomic locations of the exons. But this method isn't very scalable.
Does anyone know an easy, scalable way of retrieving the genomic locations of exon for a certain gene?

Not all exons are necessarily included in the CDS. On top of that, the first exon often is not translated completely, and only a part of the exon is part of the CDS.

Ah, sure, you are right. Somehow my brain sometimes switches the words around. Never mind me.

Perhaps you could use this (assuming hg38) as a lookup table:
$ wget -qO- https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_39/gencode.v39.annotation.gff3.gz \
    | gunzip -c \
    | grep -wF exon \
    | convert2bed --input=gff --attribute-key="gene_name" - \
    > gencode.v39.exons.bed

If you prefer to search for an Ensembl accession ID, you can remove the --attribute-key="gene_name" option. Otherwise, this puts the HGNC symbol into the ID column.
The preceding step takes about 30 seconds. Once you have made this lookup file, you can search it for exons of any gene of interest. Using MYC as an example:
$ awk -v FS="\t" -v OFS="\t" '($5=="MYC")' gencode.v39.exons.bed > answer.bed

You can then further filter that by transcript or other attributes, etc.

I'd use biomart on ensembl.org


text:
How does DESeq2 "collapseReplicates()" function work on read counts data?
Comparing read counts from an RNA-seq experiment for two select genes before and after using DESeq2's collapseReplicates() and plotCounts() functions yields interesting results:
Before collapseReplicates() and plotCounts():

Geneid
foo1.1
foo1.2
foo2.1
foo2.2
bar1.1
bar1.2
bar2.1
bar2.2
baz1.1
baz1.2
baz2.1
baz2.2
baz3.1
baz3.2

WASH7P
6
5
0
2
1
1
8
5
0
0
0
0
0
0

SOX3
1880
1861
1950
2055
1189
1181
2415
2482
3887
3810
1851
1738
3217
3406

After collapseReplicates() and plotCounts():

Geneid
foo1
foo2
bar1
bar2
baz1
baz2
baz3

WASH7P
9.877191
2.279384
3.478891
11.613875
0.500000
0.500000
0.500000

SOX3
3189.598
3563.717
3530.486
4187.011
6473.122
7991.460
5390.721

Note: In the above comparison, there are two (2) technical replicates (1.1 + 1.2, 2.1 + 2.2, etc.) for each biological replicate (foo1, foo2, etc.) for each of three (3) conditions (foo, bar, and baz).

Comparing the tables above, it appears as though DESeq2 is NOT taking the average or sum of columns being collapsed.
It is also curious - and mildly concerning - that some very low expression (i.e., 0.5 counts) is reported for genes in the matrix of collapsed replicates when, in the original count matrix, zero (0) reads were counted as 'mapped to that gene'.
So, how does collapseReplicates() "combine counts into single columns of the count matrix" as is described in the DESeq2 vignette?

Here is the code to collapse replicates and retrieve the number of read counts for a specific gene (e.g., WASH7P) in a dataframe to be used in a count plot, getting the "after" results shown above:
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ condition)
ddsColl <- collapseReplicates(dds, dds$sample, dds$run)
keep <- rowSums(counts(ddsColl)) >= 10
dds <- ddsColl[keep,]
dds$condition <- relevel(dds$condition, ref = "baz")
dds <- DESeq(dds)
countsdf <- plotCounts(dds, gene="WASH7P", intgroup="condition", returnData=TRUE)

It’s the sum. Please show reproducible code. https://github.com/mikelove/DESeq2/blob/master/R/helper.R#L186

Hmmm... alright. I have updated the question with the code (at the bottom). Let me know if you need more, like the input count matrix (cts) or annotation file (coldata).

I just realized reading the "Plot counts" section of the DESeq2 vignette that the function plotCounts "normalizes counts by the estimated size factors (or normalization factors if these were used) and adds a pseudocount of 1/2 to allow for log scale plotting", which would explain the "0.5" I observe for the genes with "0" reads mapped to those genes. This is good to know, but I now wonder if this data can be used to accurately illustrate the expression of these genes in a plot?

The DESeq2 function collapseReplicates sums the counts for the technical replicates. Here is the code reference:
github.com/mikelove/DESeq2/blob/master/R/helper.R#L186
OPs actual confusion was with the DESeq2 function plotCounts which by default normalizes count data and adds a pseudocount of 0.5 for plotting on log2 scale.


text:
Can 5' end 10x be used instead of 3' end
I have PBMCs of some patients who have recovered from cancer, and I want to look at their immune landscapes.
I got confused when looking at available kits: can I use 5' end kits from 10x instead of 3' end based?
I tried to search for information on this, but I am not clear if 5' can substitute for 3' or if they are just too different.

I'm not entirely clear on your question, but if it is about the 10X 3' vs. 5' solutions then either can be used to sequence PBMCs. The principle difference is that 3' captures the 3' end of the transcript and 5' captures the 5' end. The differences are more important for multi-modal single cell data, e.g. only the 5' can be combined with VDJ sequencing and for CITE-seq, different antibody panels are compatible with each tech. There are also some more subtle differences and technical nuances.

It depends on the task. If you want to recover TCR and BCR repertoires along with gene expression - 5' is a must. Also, for an overview repertoire analysis you dont even have to perform target enrichment and use 5' RNA-seq daata to directly recover the repertoire.
Its easy to do with MiXCR software in one line command:
mixcr analyze 10x-5gex-full-length \
 --species hsa \
 sample_R1.fastq.gz \
 sample_R2.fastq.gz \
 sample_result

Though MiXCR currently only supports repertoire analysis it will also support gene expression soon.
More on that here:
https://docs.milaboratories.com/mixcr/reference/overview-built-in-presets/#10xgenomics
Disclaimer: I work for the company that develops MiXCR.

Just to mention that it is possible Mark Izaelson is employed by MiLaboratories. However, it is great that the post was edited to describe the algorithm.


text:
#1 Parsing trouble java.lang.NumberFormatException in GSEA 4.2.1
I am getting the following error message when I tried to load .TXT data file to GSEA 4.2.1.  I tried different way to change the file name as .GCT but not working.  Any suggestion? Many thanks.
Full Error Message:
There were errors: ERROR(S) #:1
Parsing trouble
java.lang.NumberFormatException: ...

---- Stack Trace ----
# of exceptions: 1
------For input string: "description"------
java.lang.NumberFormatException: For input string: "description"
    at java.base/jdk.internal.math.FloatingDecimal.readJavaFormatString(Unknown Source)
    at java.base/jdk.internal.math.FloatingDecimal.parseFloat(Unknown Source)
    at java.base/java.lang.Float.parseFloat(Unknown Source)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.AbstractParser.parseStringToFloat(AbstractParser.java:266)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.AbstractParser.parseFieldsIntoFloatArray(AbstractParser.java:371)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.AbstractParser.parseTextMatrixToDataset(AbstractParser.java:293)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.TxtDatasetParser.parse(TxtDatasetParser.java:109)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.readDatasetTXT(ParserFactory.java:201)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:760)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:737)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserWorker.doInBackground(ParserWorker.java:52)
    at java.desktop/javax.swing.SwingWorker$1.call(Unknown Source)
    at java.base/java.util.concurrent.FutureTask.run(Unknown Source)
    at java.desktop/javax.swing.SwingWorker.run(Unknown Source)
    at java.base/java.util.concurrent.ThreadPoolExecutor.runWorker(Unknown Source)
    at java.base/java.util.concurrent.ThreadPoolExecutor$Worker.run(Unknown Source)
    at java.base/java.lang.Thread.run(Unknown Source)

Can you elaborate on what was not correctly formatted? I followed the instruction on their manual .txt and am getting this same error. Thank you.

Sorry for my delayed response, if you do get chance @Hong please remember to "accept" your answer. Much appreciated.

I just figured out the error. My .TXT was not correctly formated, following the instruction from their manual under .TXT format worked.


text:
In a RNA-Seq heatmap should you do Z-score standardisation before clustering the rows/columns or after?
I have made a heatmap using RPKM values from a RNA-Seq dataset using the pheatmap() function in R. I have log2-transformed the data before performing Z-score standardisation of the data. I have also clustered the rows and columns of the heatmap. I have the following code:
heatmap_data_log2 %>% pheatmap(color = colorRampPalette(c("blue2","white","red"))(100), scale = "column", cluster_cols = T, clustering_method = "ward.D2", angle_col = 45, fontsize_row = 7.5, fontsize_col = 8, border_color = NA, cutree_rows = 4)

I have seen in a video tutorial that the pheatmap() function does the Z-score scaling before doing the clustering of the rows and columns. However, in this online article, it says that "The Z-scores are computed after the clustering, so that it only affects the graphical aesthetics and the color visualization is improved".
As these two sources are giving contradictory information, I was wondering which is better. Should you do Z-score scaling of the gene expression values before doing the clustering or after? Any advice is appreciated.

Looking at the source for pheatmap, there is a function called scale_mat that is used to preprocess and normalize the input matrix, depending on the value of scale, which specifies one of either none, row, or column normalization options. Separate and similar functions are used for color value scaling downstream.
I have not watched the video or read the linked article, so I'm not sure this is contradictory, but it may be just small confusion about how normalization is applied (and to what).
Normalization is usually done before clustering to focus cluster construction on signal rows or columns that are not some multiple or scale factor of another — and in referring back to the source code, this is what pheatmap does, as well.


text:
How do I set my path for the SRA toolkit under macOS 12.1?
I would like to add the SRA toolkit to my path under macOS 12.1.
After downloading the SRA directory, decompressing, and adding to the path, the directory does not seem to be available e.g. I'm attempting to add it to my path with
export path=(/Users/ray/Documents/School/bioinformatics/sratoolkit.2.11.2-mac64/bin)

and it does not work. For example if I type
 `ray@dynamic-xx-xxx-xxx-xx ~ % which prefetch

prefetch not found` macOS gives me prefetch not found.  The result of my path is
`ray@dynamic-10-125-206-18 ~ % echo $PATH
/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin`

Is there a way to add it to my path so it works?

Firstly, you are trying to set a path into the Documents directory? Thats like for documents.

Anyway there could be a permission issue. The quickest way is simply set a path within .bashrc or .profile whatever into ~/bin
and then use
ln -si /Users/ray/Documents/School/bioinformatics/sratoolkit.2.11.2-mac64/bin/prefetch ~/bin

... for each compiled program.
$ prefetch

should then work.
This creates a symbolic link in ~/bin so when you type the program in your commandline prompt it will fire, without wiring it via a custom bin. If there is a permission issue a permission error will be quickly flagged.

The best way is to use Anaconda because its within the Anaconda collection and simplifies everything on three counts, viz:

in answer to your question all configurations are automated
it is easy to automate all updates
it offers a virtual environment where code can't mess with your system
whenever you want to download other packages its the same approach

Thus all you have to do is download and install Anaconda or miniconda. From this point everything is done using conda, either:
conda install -c bioconda sra-tools
conda install -c bioconda/label/cf201901 sra-tools

If one fails the other one will work. Then simply,
conda activate

SRA is now fully configured and ready to go. To get out of the conda virtual environment,
conda deactivate

If you use Anaconda alot it is better to make a specific environment for SRA and not use the core environment, but thats another question.
Anaconda installs in two different locations depending on whether you install the GUI version, which makes maintenance and environments a bit easier. In the later case I think its /opt (but forget), otherwise the commandline version is ~/anadonda3.

I see your using a virtual environment via homebrew. Anaconda is separate, its the same idea, but its better.

Fixed it by installing anaconda. It does look like sra-toolkit is not compiled for Apple Silicon yet so Apple x86_64 miniconda is the only architecture that sra-toolkit will work with on MacOS.

Ah right, the M1 chip compatibiility issue. Good to hear that Anaconda has got around the issue.


text:
what are the nodes and edges values in gene regulatory networks?
I am trying to find out how one can using gene expression data can infer gene regulatory network applying graph theory concepts. But I could not find a proper reference that

explain how one can get the adjacency matrix for a gene regulatory network? like what are nodes what are the edges? which kinds of parameters one should infer?

how one can use gene expression data for this modeling?

Any reference or help regarding these questions is highly appreciated if even it is going to be a keyword.

Hi @Raz, this does not seem like a specific question, and needs more information to be properly answered. Could you please add in information about the specific situation or problem that you are trying to solve?

This is an incomplete answer since it's not my field of expertise, but it sounds like Weighted correlation network analysis is what you are referring to.
There is a paper which describes how it can be applied to gene expression data.
Langfelder, P., Horvath, S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9, 559 (2008). https://doi.org/10.1186/1471-2105-9-559

There are many different kinds of "gene regulatory networks", so how to model them with gene expression data depends on what you're trying to model and what type of regulatory mechanisms you want to study.
Pick a random gene in the human genome and use a resource database like StringDB.
Let's take the TET1 protein for example.
Some edges in that graph represent co-expression that you'd be able to find from protein expression data (e.g. NANOG, DNMT3B, and more are linked to TET1 this way).
This type of graph tries to model how different proteins interact with each other inside a cell.
Alternatively, you can look at 3D chromatin interactions from Hi-C data, for example.
Nodes in that graph could be regulatory elements like enhancers and promoters, and the edges could be how frequently these elements interact with each other.
This type of graph tries to model how DNA elements interact with each other to control expression of nearby genes.
This is a different purpose to the one above, so it needs to be modelled and interpreted differently.
Each of these are their own type of "gene regulatory network" which require their own methods for constructing adjacency matrices.
For example, you could construct an adjacency matrix by looking at the correlation of expression for each pair of proteins (or transcripts, or genes) across a large set of samples.
Or in the promoter-enhancer interaction case, a Hi-C contact matrix is itself a type of adjacency matrix.
Different graph models (and thus different adjacency matrices) will be useful for different purposes.
So it might be worth trying to refine what you mean by a "gene regulatory network", what aspect of biology you're trying to study using this network, and how different types of gene expression data may be useful to construct these networks.
John Quackenbush's and Andrea Califano's labs have produced a lot of papers modelling different types of gene regulatory networks, like The Network Zoo and ARACNe.
These mostly use gene expression data from RNA-seq measurements, but you can read up on some of these to get some ideas on how these types of models may work.


text:
How to input data and metadata from NCBI for RNA-Seq analysis in R
So, I have downloaded some .txt and .csv files containing the raw count data (I am supposed to find out differentially expressed genes and see which genes and up- and down-regulated). I have also downloaded both the .family.soft and the series matrix files as they mentioned they contained metadata - however I'm facing several errors while trying to import them.
And all the tutorials I am trying to follow seem to have premade metadata files and there is no explanation on how this metadata is made.
Please help me in trying to understand what I am doing wrong and how to proceed. I will try to give as much information as I can to the best of my abilities, as I have barely any experience in this field.

You may want to add to the post the first ten lines (and if the files are delimited with >10 columns, restrict to just the first 10 columns as well) of each file so people here can see what these files contain.

A metadata file isn't that large or complicated.  If you want to do DE, you must have two treatments?  Or two timepoints?  Two or more something, and you are comparing samples of one type to another?
It's just a data table where the rownames are the samples, and then a column for, say, treatment.  If you have some samples treated, some not, some at day 0, and some at day 5, then you have a column for Day and another column for Treatment.  (Depending on what you want to compare to what, you might want another column of Day_Treatment)

I don't expect the soft or matrix files will be useful for you, although it might be possible to pull metadata out of those files if absolutely necessary. I usually only see those for microarray data, so it might be worth double-checking that you do actually have RNASeq data.
I've found the best "quick start" explanation on how to carry out differential expression analysis to be the one for DESeq2. As input, this requires three things:

A count matrix, with samples in columns, genes in rows, and matrix entries being the raw read counts of reads mapped to genes for each sample
A metadata data frame, with rows arranged exactly as in the columns of the count matrix, and fields describing features of each sample
A statistical model / equation that describes the relationship between samples

I have found the equation to be the most complicated to understand, and unfortunately it's very situation dependent. In its simplest form, setting up a DESeq2 data structure would look like this:
dds <- DESeqDataSetFromMatrix(countData = count.mat,
                              colData = metadata.df,
                              design= ~ Condition)

This would be for a situation where all samples were sequenced in the same sequencing run, and only a single characteristic were being compared (e.g. with vs without treatment). The metadata data frame might look something like this:
metadata.df <- data.frame(row.names = colnames(count.mat),
                          Condition = c(rep("treated", 6),
                                        rep("untreated", 6))

e.g.:
          Condition
Sample01    treated
Sample02    treated
Sample03    treated
Sample04    treated
Sample05    treated
Sample06    treated
Sample07  untreated
Sample08  untreated
Sample09  untreated
Sample010 untreated
Sample011 untreated
Sample012 untreated


text:
Acquiring & Processing single cell data (not 10x)
I have a question related to processing of Single cell RNA-seq data originating from non 10x platform. I have worked with data that originated from 10x platform and I parse it through cell ranger pipeline and then put it through Seurat.
As an example the dataset in this GSE156644 repository are not being processed through Cell ranger pipeline and I get an error that indicates "single-end read found" another example is of this GSE115235, in the later I have bunch of files that look like bulk RNA-seq.
If you check the GSE115235 it has bunch of files that are not cellranger compatible. However, they are from single cell experiment. I find it strange that NCBI accepts such submission and the reviewers don't even pay attention to what has been submitted.
I guess, that these are SE RNA seq data that can be aligned using STAR or Bowtie and then merged in some way to get per cell per sample expression.
Along similar lines, how do I deal with 2 fastq files for a scRNA seq run? example https://www.ebi.ac.uk/ena/browser/view/PRJNA602526
This link has scRNA seq files that look paired end (R1 & R2 but NO I1). I am bit confused as cell ranger count option requires the folder to have 3 files (I1,R1,R2). How can I run cell ranger with R1 & R2? Or do I need to use another aligner?
I will send an email to sra and also to the authors, but any other help to navigate through these datasets and process them is deeply appreciated.
Thank you.

Cell ranger does not require I1 and it was a bad idea for them to ever support that. You can directly run cellranger count on demultuplexed fastq files. 2. Salmon alevin or STAR-solo would be examples of common alternatives. Salmon alevin is in fact much faster and likely better.

Modified from @devon-ryan's comments:
The couple of samples I looked at have the sequence from read 1 embedded in their header. I1 isn't needed, they're already demultiplexed. You could feed those into STAR or similar tools, but processing them will be more of a pain since most tools expect the original unprocessed version of what they uploaded.
Modified from atpoint's comments:
Single-cell data is not really useful as "bulk" as the notable amount of PCR during the library prep creates quite a bias that without the UMI information cannot be corrected. I'd email the authors and ask for the R1 files. Single-cell uploads at NCBI tend to be a mess, that's unfortunately how it is. Without R1 data are almost useless.
Those file format issues are strange, but it is what it is. Reviewers do not pay attention to method sections, choice of tools and data submission most of the time in my experience, especially not when the paper is biology-focused rather than a computational-focused one. Even more, you often do not find the suitable metadata for a proper analysis such as information about batches etc. The only choice you have is emailing the authors as noted above.
If by "two files" you mean that the same library (=sample) was sequenced over multiple lanes (very common to do that since required depth is often larger than what a single lane can provide) then merge the fastq files, so all R1s into one and all R2 into one file.


text:
How to improve a genome assembly using Dovetail and PacBio assembly?
I have more of a conceptual question. I have two genome assemblies from the same plant, one from Dovetail technology (~998 Gb) and another is PacBio HiFi assembly (~1.1 Gb). The Dovetail assembly is more contiguous but has lower base quality whereas the PacBio has higher base quality but is more fragmented. It is a diploid organism with relatively high heterozygosity.
I tried to use RagTag patch to improve the assembly. But no improvement was obtained.
Is there a way to use both assemblies and produce a hybrid assembly with high contiguity and base quality?

You can't do assembly with Dovetail alone. You have to have WGS reads in the first place. Is the Dovetail assembly derived from HiFi assembly?

Answer from @maximilian-press, converted from comment:
I am confused when you talk about base quality, because those technologies are largely about ordering and orienting existing contigs, not assembling contigs in the first place. If you have Hi-C reads, you can read them into juicebox and actually visualize the junctions to see if the assembly makes sense: aidenlab.org/juicebox. It would also be helpful to see the actual numbers, not just "more" or "lower".

But do not you think that it is possible to use the Dovetail assembly to scaffold the PacBio HiFi assembly with more contiguity?

@AnikDutta I think that the specific purpose of "Dovetail" (which I assume means Hi-C and/or Chicago) is the scaffolding of genomes. but on its own it is not much good for assembly. A "Dovetail" assembly is probably a set of contigs created by non-Dovetail technology (like PB HiFi), which have been scaffolded using an extra layer of "Dovetail" data. To be able to answer this question we would need to know what that non-Dovetail technology is that underlies the final assembly.

@AnikDutta For example: are the PB HiFi contigs the basis of the "Dovetail" assembly? If so, then the "Dovetail" assembly probably already represents an effort to scaffold the HiFi assembly.

@MaximilianPress Thanks a lot for the reply. No, the PB HiFi contigs are not the basis of the Dovetail assembly. As far as I know, Illumina reads were used for Dovetailing. The PB was generated independently. I ran the BUSCO and found out that the PB HiFi gives 98% BUSCO score whereas the Dovetail gave 95%. What do you think of using the Hi-C data and HiFi data to create haplotype resolved assembly? Will it make sense here considering it is a highly heterozygous diploid species?

@AnikDutta If the DT assembly is Illumina-based, I would strongly recommend the workflow of HiFi + Hi-C assembly. hifiasm natively takes in both of these data types: https://github.com/chhylp123/hifiasm, it should be pretty easy to run as long as you have the compute resources. You can polish the assembly post hoc using juicer/juicebox, after mapping the Hi-C reads onto the resulting assembly, as suggested above.

Yes, that was what seems the most appropriate wayforward. Thank you very much.


text:
Is the algorithm of the STAR RNA Seq-Aligner similar to the Knuth-Morris-Pratt string matching algorithm?
I don't know if I'm in the right SE, if not let me know and I'll delete it.
I am reading the original publication of the STAR RNA-Sequence aligner.
Given the quoted text below I wonder if the algorithm for the seed search is actually an implementation of th Knut-Morris-Pratt algorithm?

The central idea of the STAR seed finding phase is the sequential search for a Maximal Mappable Prefix (MMP). MMP is similar to the Maximal Exact (Unique) Match concept used by the large-scale genome alignment tools Mummer (Delcher et al., 1999, 2002; Kurtz et al.) and MAUVE (Darling et al., 2004, 2010). Given a read sequence R, read location i and a reference genome sequence G, the MMP(R,i,G) is defined as the longest substring (Ri, Ri+1, … , Ri+MML−1) that matches exactly one or more substrings of G, where MML is the maximum mappable length. We will explain this concept using a simple example of a read that contains a single splice junction and no mismatches (Fig. 1a). In the first step, the algorithm finds the MMP starting from the first base of the read. Because the read in this example comprises a splice junction, it cannot be mapped contiguously to the genome, and thus the first seed will be mapped to a donor splice site. Next, the MMP search is repeated for the unmapped portion of the read, which, in this case, will be mapped to an acceptor splice site. Note that this sequential application of MMP search only to the unmapped portions of the read makes the STAR algorithm extremely fast and distinguishes it from Mummer and MAUVE, which find all possible Maximal Exact Matches.

It is not cited in the references but it might be considered 'common knowledge' by now which is why the citation is omitted...

No, STAR isn't using the KMP algorithm or a modification of it.  The KMP algorithm is an online exact pattern matching algorithm.  It does (linear time) pre-processing on the query and then finds all exact occurrences of the query in the string (as you may be alluding, you can also modify KMP to report maximal exactly matching prefixes of the query).  This process takes time of order N+M where N is the length of the reference and M the length of the query.
On the other hand, STAR pre-computes a suffix array on the reference, and then uses this to answer exact prefix match queries across many different queries.  This is a very different approach.  Building the suffix array takes time linear in the reference, and storing the suffix array requires space linear in the reference (though the constant factor is pretty big).  However, it allows enumerating all occurrences of a query Q in the reference R in time $O(k + lg(\left|R\right|)+ \left|Q\right|)$ where $k$ is the number of occurences of Q in R.  This holds not just for the whole pattern, but for any maximal prefix as well.  Also, in practice, STAR does this much faster by building a lookup table that maps small prefixes (say, all length 12 prefixes) to the contiguous suffix array interval where they occur.  In practice, this gives very fast lookup of queries.  Again, a key difference between this kind of approach and KMP is that the suffix array pre-processes the reference and therefore builds an index that works for many different queries (the sequencing reads), while KMP pre-processes the query to find all of its occurrences on the reference.  Another key difference is that, while the KMP algorithm is always linear in the length of the reference to return all matches in the worst case, the suffix array can be considerably faster.
In the original STAR publication, this is mentioned in the supplementary material in section 1.2 "Pre-indexing of suffix arrays".  This idea (along with others) is also covered in quite some depth in this paper.  Basically, the look up table has been a folk technique for a long time, but it's not clear there is an "original" source for it.

a sequel in case you are interested in answering it: https://bioinformatics.stackexchange.com/questions/18388/how-should-i-imagin-the-pre-indexing-of-suffix-arrays-in-the-star-rna-seq-aligne


text:
How to interpret the Copy Number from pseudocount depmap portal file?
This question has also been asked on Biostars
How do I convert the Copy Number (CN) file into useful information such as HOMDEL and AMP?
This is an already normalized file with pseudo count + 1. In this file there are values from 0 to 1 for each cancer cell line, in wide format.
I have already looked at the forum post about the file, but it's not helpful for me.
This is what the file looks like:
,DDX11L1 (100287102),WASH7P (653635),FAM138A (645520),OR4F5 (79501)
ACH-000621,1.0166833616923812,1.0166833616923812,1.0166833616923812,1.0166833616923812
ACH-002512,2.9659871795507433,2.9659871795507433,2.9659871795507433,2.9659871795507433
ACH-001188,0.9881267634909964,0.9881267634909964,0.9881267634909964,0.9881267634909964
ACH-000471,0.8053837915961594,0.8053837915961594,0.8053837915961594,0.8053837915961594

Thank you for your help

I'm not familiar with this file structure, but have you looked at the forum post about relative copy number linked to in the FAQ?

Since we do not have matched normals, the output is a “copy ratio” or relative copy number. It is relative to the rest of the genome for that cell line. E.g. if the cell line is tetraploid we would not be able to see it from the relative copy number. These values are reported as log2(relative CN + 1) in the portal.

Yes I checked but I still have the same doubts


text:
How to extract and visualize sequence logos from CNN kernels?
I created a basic CNN architecture using Tensorflow to classify transcription factor binding sites. My aim is to somehow extract and visualize sequence logos from the convolutional kernels. The model looks like this:
model = tf.keras.models.Sequential()
model.add(tf.keras.layers.Conv2D(128, (1, 24), padding='same', input_shape=(1, 101, 4), activation='relu'))
model.add(tf.keras.layers.GlobalMaxPool2D())
# ...

I am using the supplement dataset from this publication
What I tried so far is that I have trained the model and created a new model out of it which has only an input layer and an output layer that contains the weights of the original convolution kernels:
test_model = tf.keras.models.Model(inputs=model.inputs, outputs=model.layers[0].output)
feature_maps = test_model.predict(x)

I have ran a prediction on a sequence and plotted out a convolution kernel:

At this point I got a bit stuck as I am not sure how can I move forward with this in the direction of visualizing the sequence logos.

What you are missing is the accuracy of the prediction, this is the next step and must be >80%, preferably >85% and >90% is good in biological systems. If its lower than 80% you'll need to adjust the training. Following this there are issues with false positives and false negatives.
The method of using seems fine for mapping sequence motifs to transcription factors, which is my interpretation of your question. There must have been vectorisation of the sequence data and this is not trivial.
To answer your question, which you might be asking, you cannot produce an automated heatmap of sequence motif in relation to the transcriptor factor binding. You simply input a sequence and the associated transcription factors are revealed. Unless you vectorised your transcription factors somehow and associated this with a sequence motif.
To produce a heatmap, which again is what I suspect you are requesting, you do this via iteration. So you need a bit of coding. Thus you take all relevant sequence motifs - particularly those that have not been trained - and input them into deep learning algorithm and map which motifs associate with which binding factors and assign a numerical score to given groups of sequence motifs. Any graphics program can then be used to convert this to a heat map. Deciding 'taxonomic groups' for the sequence motifs is subject to your discretion. If you've given the transcription factors and trained this onto the sequence motifs, then classification of transcription factors I suspect would be much easier.
Heatmaps are produced in other machine learning approaches, i.e. which exclude deep learning, but your next step is percentage accuracy regardless of approach.


text:
PSSM Matrix in PSI-BLAST
After running PSI-BLAST one obtains a profile matrix with this header:
 A   R   N   D   C   Q   E   G   H   I   L   K   M   F   P   S   T   W   Y   V

How does PSI-BLAST handles X and B "residues" in the fasta file then?

X residue being an amino acid that can't be assigned due to the presence of nucleotide ambiguities which result in the codon representing more than one amino acid.
B (beta) residues, with the exception of B-alanine are singly generated in vitro and being explored for their putative antibacterial properties.
PSSM is the position specific score matrix

It's a blast thing in general. Degenerate bases are ignored in DNA queries. Pyr, asx & glx in protein also. I don't think they count as a gap, but instead these IUPAC letters and any random letter (=non-space/non-punctuation) character (even Unicode IIRC) become a masked character —they will be in grey if so: it would straightforward to test by simply doing a pairwise alignments with one sequence with asparax (B) and glutamax (Z) replacing the respective amide/acids. There should be written somewhere how masked characters behave during scoring.


text:
How should I imagine the pre-indexing of suffix arrays in the STAR RNA Seq-aligner?
I'm using STAR for an internship (fusion-genes in ALL cancer) and I'd like to understand what exactly it is doing. The people at the internship don't know either how it works, given their biology backgrounds. They told me "its enough to know that it works". That was not an satisfactory answer for me which is why I am fighting my way trough the suplementary materials trying to wrap my head around the algo. My background is bio as well but I've spent reading a book about algos and algos in BI this year.
I am trying to understand the section Pre-indexing of suffix arrays from supplementary materials of the Dobin et al. STAR RNA seq-aligner Publication.

While suffix array search is theoretically fast owing to its binary nature, in practice it may suffer from non-locality resulting in persistent cache misses which deteriorate the performance. To alleviate this problem we developed a pre-indexing strategy. After the SA is generated, we find the locations of all possible $L-mers$ in the SA, $L<=L_{max}$, where L_{max} is user defined and is typically 12-15. Since the nucleotide alphabet contains only four letters, there are NL=22L different L-mers for which the SA locations have to be stored.

The authors also state:

Suffix Array (SA) of the whole genome is utilized to find the Maximum Mappable Prefixes (MMP).

I am confused by the use of the term $L-mer$. When googleing I only keep finding $k-mers$ in combination with suffix arrays. I did however find this GitHub page mentioning $L-mers$ in the following way:

This program computes either the suffix array or the L-mer array of a string. The suffix array gives the lexicographical order of all the suffixes of a string. The L-mer array gives the lexicographical order of all the L-mers in the string (sorts the suffixes with respect to only the first L characters).

→ Can I assume that $L-mer$ was used instead of $k-mer$ to stress the lexicographical ordering? But as far as I understood all suffix arrays contain the suffixes' index based on its lexicographic ordering as illustrated by wikipedia's page on suffix arrays:

Do I have to imagine the entire genome being a single string that is converted into one single suffix array? Could you provide sources fit for a beginner with biology background and very little idea about algorithms?

According to BI SE this falls into the category ".... However, if your motivation is “I would like others to explain ______ to me”, then you are probably OK."
From https://bioinformatics.stackexchange.com/help/dont-ask
Which is why I was hoping to find some help here.

Thanks for the clarification. It's fine to ask these types of specific questions, but context helps a lot in providing useful answers. You've mentioned "a beginner with biology background and very little idea about algorithms", which is inconsistent with the level of algorithmic detail you've provided. That's why I asked for clarification.

@gringer Thanks, I see (-:  trust me, algorithms are still a big challenge and mystery to me.

You can take a look at slide 40 of the lectures slides I use when I teach suffix arrays (reproduced here)

An L-mer is a K-mer — STAR builds an explicit lookup table from all length 14 (by default) strings that maps the length 14 string to interval of the suffix array that contains all suffixes that begin with this length 14 prefix.  This means that, to find a MMP, instead of having to do a binary search on the entire suffix array, you only have to do a binary search on the sub-interval that begins with the length 14 prefix that your query string starts with.  If all substrings suffixes began with length 14 prefixes with equal frequency, then this would reduce the size of the interval you have to look at by a factor of $268435456 = 4^{14}$.  Of course, in biological sequence, all prefixes are not equally likely, so the speedup is more variable in practice.  Nonetheless, it makes lookup in the suffix array much, much faster.
Of course, in addition to the core algorithm, STAR contains many heuristics not all of which are likely described in the paper (especially since the tool has been and continues to be under active development and improvement since the paper's publication).  Nonetheless, this is the idea behind the suffix array pre-indexing.  Since STAR still has the whole suffix array, if there is no length 14 exact prefix match, it is still possible to find a smaller match using a more standard binary search in the suffix array.  However, I do not know (and am not sure where it is documented) exactly when/if STAR will fall back to a more naive binary search if a prefix L-mer is not found.  If you are truly curious about this, I'd recommend asking Alex Dobin (the STAR author and maintainer) over on the STAR GitHub repository.  He is generally very responsive to user queries there.


text:
How to consume non-CRAN packages in custom package?
In a custom package which imports a package from Bioconductor, what is the right DESCRIPTION file?
I added in my DESCRIPTION a line with Remote like so
Remotes:
    bioc::3.14/dada2

as suggested by https://cran.r-project.org/web/packages/remotes/index.html vignette, but it still failed the R check() with error "Packages required not available". I am very confused about what is the right way to solve this problem, since multiple sources propose a different solution:

https://cran.r-project.org/doc/manuals/r-release/R-exts.html#Package-Dependencies mention Additional_repository key in the DESCRIPTION

mention Remotes, suggesting one add a Remotes entry to DESCRIPTION with a specific syntax for Bioconductor dependencies

The following post suggest using blockViews R package development: How does one automatically install Bioconductor packages upon package installation?

What is the appropriate way to declare dependencies on packages on Bioconductor? The remotes support helped me satisfy a dependency on a package from Github, but not for bioconductor

As I answered on the linked question, the right way to declare packages from repositories recognized by R is to just declare normally on the Imports, Depends, Enhances and LinkingTo.
Bioconductor is an accepted repository so it doesn't need to be added on Additional_repository field.
BiocViews is a required field from Bioconductor. As CRAN accepts other fields on the DESCRIPTION it could be used by CRAN packages.
You don't mention it but if you found this problem and it wasn't on your computer were you are developing the package it means you are using a continuous integration test. The problem with CI like Github actions, travis, appveyor or other is that dependencies must be installed before (as they should on your computer).
To install these dependencies before that depends on the CI you use and the strategy used there (install from ppa, use RSPM, use bspm::install, ...) . However, I recommend to do it via BiocManager::install because this way you can install all packages from both CRAN and Bioconductor and their dependencies.
Even if the package is on CRAN, github or just Bioconductor with BiocManager::install you'll be able to install all the dependencies and continue to check your package normally on the CI.

Additional_repositories is for other non standard repositories like a personal one or other a specific of a project. It is used by only ~100  packages of the 18765 packages on CRAN. But there are many more dependencies from CRAN to Bioconductor (and vice versa) .

That's a surprising behavior.  I need to install BiocManager manager and all the Bioconductor dependencies to succeed in check(). This is not really very user friendly. Is there a way to do that with a single command in a fresh R environment?

This is very user friendly, as dependency installation (of the right Bioconductor release) is taken care for you. Otherwise you need to manually chase down each dependency from each repository on the right version. In a single line: R -e install.packages("BiocManager"); BiocManager::install("my_package", dependencies = c("Depends", "Imports", "LinkingTo", "Suggests", "Enhances")) This will install my_package and all the dependencies.

A note about check: if you are using check you are not a user, you are a developer of R extensions. To be able to provide high quality packages to users R CMD check my_package is used, but to check the package it needs all the dependencies (otherwise it cannot fully check the package). If you don't want to check your package you can simply build the package R CMD build my_package and share the tar.gz file. This will provide the package to users without any check (it could even be that your package is not installable). Hope this helps!

Thank you. Is there a way to install all dependencies from DESCRIPTION you are suggesting without running an install of my_package? Something that I need to run before R CMD check on a new renv. This is very user friendly, as dependency installation (of the right Bioconductor release) is taken care for you. Otherwise you need to manually chase down each dependency from each repository on the right version. In a single line: R -e install.packages("BiocManager"); BiocManager::install_deps(dependencies = c("Depends", "Imports", "LinkingTo", "Suggests", "Enhances")) or something like that?

If you are using Github actions I recommend to check out biocthis as it helps setting this. You could also use remotes::install_deps(repos = BiocManager::repositories()) which would only install dependencies from the right repositories but not the package being installed.

Thanks, very useful ! Can you share with me the links to the relevant documentation?

yes, you can find them  via ?remotes::install_deps, ?BiocManager::repositories on the terminal. And here are links to their respective manuals: remotes, BiocManager


text:
Polygenic Risk Scores
I have run a GWAS and used PRSice (https://www.prsice.info/) to calculate polygenic risk scores on the data collected from the GWAS. I need to check whether alcohol dependence can predict gambling disorder, through the PRS. However I am having trouble with how exactly to do that. I have a GWAS on alcohol disorder from the one population and a GWAS on gambling and alcohol disorder from the other population, then I calculated PRS for the three GWAS datasets, but how do I use the one to predict in the other? Through the literature it was suggested to use linear mixed models, which I have done but I still feel as if I am doing it wrong as I only did the model on the one dataset and not the two combined.

Before doing anything else, you should split your dataset into training and test groups, and [this is important] don't touch the test dataset until you're fully confident that your approach is correct: no post-test modification.
A split that's considered acceptable in Computer Science is a 10%/90% split, i.e. use 10% of the data for training / algorithm generation and the remaining 90% for testing. This will give you good confidence that the algorithm works.
I personally think it's a good idea to do population sub-sampling (i.e. further subsample that 10%) to determine an informative marker set prior to any risk calculations. What you should be looking for are markers that are consistently linked with the trait of interest (i.e. in your case alcohol dependence) in multiple sub-sampled populations. I go into a lot of detail about how to do this here [bearing in mind it's an old pre-print, and doesn't do the 10/90 split I've recommended above].
After you're happy with the marker set, calculate the risk scores only using that set (hopefully PRSice can do this), then carry that risk score [only] through to test against the dependent trait of interest in the training dataset (i.e. in your case gambling disorder). If your dependent trait is a binary trait, this could be done using ROC analysis (as in my preprint) to calculate area under the curve. Or use linear mixed models, whatever works for you.
Make any changes you think you need to improve the test result (e.g. changing the composition of the marker set).
After convincing yourself that the approach is correct, do the same risk score calculation (using the same set of markers) on the test set, find out how well it predicts the dependent trait, and report on the results. Don't change anything.


text:
Assign multiple taxids to a sequence when constructing a local BLAST database
I recently had a script fail due to poor handling of BLAST output. The BLAST -outfmt staxids field usually returns a single taxid, but occasionally it returns two or more taxids separated by a semicolon, such as 556514;701533. Fixing the script to handle this should be fairly straightforward. But the original failure happened when I was searching a local copy of the NCBI nt database, and I'd like to create a test case with a small query and small database that reproduces this issue more rapidly.
I am familiar with the -taxid and -taxid_map flags of the makeblastdb command, but I have been unable to reproduce this behavior with BLAST; that is, make BLAST return multiple taxids for a sequence.

If I provide 556514;701533 or 556514,701533 as a sequence's taxid with -taxid or -taxid_map, the makeblastdb command fails citing an integer conversion error.
If I provide 556514 701533 as the taxid in the taxid_map file, the database will be indexed successfully but BLAST will return the invalid value 0 for staxids.
If I duplicate each entry in the -taxid_map file so that there is an entry for each taxid (see below), only a single taxid is reported by BLAST.
seq1 556514
seq1 701533
seq2 556514
seq2 701533
seq3 556514
seq3 701533
...

When constructing a custom local BLAST database, how do I associated a database sequence with more than one taxid?

UPDATE: The offending entry (or one such entry) appears to be a duplicate or merged entry in the nt database. When I use blastdbcmd to retrieve the subject sequence by accession, nt returns a sequence with two deflines on single line, like so.
>accession1 plain text description >accession2 other description
ACGT...

Using either of the two accessions, blastdbcmd will return the exact same entry, and the two staxids shown above correspond to the two distinct accessions.
So part of the issue seems to be how to construct a BLAST database such that multiple accessions can point to a single entry/subject.

Can't you just jiggle it around in pandas? That was my approach leveraging local blast, albeit for something it was probably not intended to do

@M__ I already have what I think is a solution, yes involving post-processing the BLAST tabular results in a dataframe. But I'll re-emphasize that I'm trying to construct a test case that will reproduce this behavior in BLAST.

Thanks now I understand. If a Blast expert answers I'll be enlightened too.

I inquired with the NCBI BLAST help desk today, and within an hour they provided a very helpful response.
To summarize, when redundant records in the database are collapsed, the deflines (accessions + descriptions) of the original entries are placed on a single line separated by an invisible control character (Ctrl-A or \x01). The makeblastdb program will recognize this and associate multiple accessions (from the combined defline) and taxids (from the taxid_map file) to a single entry. If you are retrieving entries from e.g. nt or nr with blastdbcmd, you can use the -ctrl_a flag to format the data with these invisible control characters—by default the > character is used.

tl;dr I can reproduce the behavior described in the original question by:

creating a FASTA file containing an entry with multiple accessions separated by the invisible \x01 character
create a taxid_map file mapping the accessions to different taxids
format the FASTA file with makeblastdb -taxid_map ...
perform a BLAST search with a query that matches the entry in question; -outfmt "6 staxids" will return a string with multiple taxids.

Thanks for answering your question. Its good to know the NCBI helpdesk is fast and knowledgeable.


text:
Different output on same command and same input file, different library version
I am running a variant calling pipeline to detect solid tumours and haematological alterations. The variant calling is done with Mutect2, whereas I am using an R script to perform some filtering on the VCF produced by Mutect
I noticed, however, that the output of the header command (VariantAnnotation library, R programming language) is not the same considering two versions (1.24.5 vs. 1.32.0).
When using v 1.32.0 I get
> header(input)
class: VCFHeader
samples(1): Panel-STHT-T210986-PJR
meta(6): MutectVersion fileformat ... GATKCommandLine contig
fixed(1): FILTER
info(16): CONTQ DP ... STR TLOD
geno(13): GT AD ... SAAF SAPP

Whereas when using v 1.24.5 I get fixed(0): instead of fixed(1): FILTER, as shown above:
> header(input)
class: VCFHeader
samples(1): Panel-STHT-T210986-PJR
meta(3): META GATKCommandLine contig
fixed(0):
info(16): CONTQ DP ... STR TLOD
geno(13): GT AD ... SAAF SAPP

I also noticed that if I continue with the filtering that the R script does (eventually writing a vcf with writeVcf), with VariantAnnotation v 1.32.0 I get an extra line in the vcf header: ##FILTER=All filters passed, which is upsetting another downstream program (the error message that I get is .vcf: Invalid FILTER line). And it appears to me that this ##FILTER=All filters passed line is actually coming from the extra line of the header command.
The easiest solution for me would be to downgrade from 1.32.0 to 1.24.5: However, I would like to understand where this extra line comes from
Session information (when using v1.32.0 of VariantAnnotation is shown below:
R version 3.6.1 (2019-07-05)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: CentOS Linux 8 (Core)

Matrix products: default
BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so

locale:
 [1] LC_CTYPE=en_GB.UTF-8          LC_NUMERIC=C
 [3] LC_TIME=de_BE.UTF-8           LC_COLLATE=en_GB.UTF-8
 [5] LC_MONETARY=de_BE.UTF-8       LC_MESSAGES=en_GB.UTF-8
 [7] LC_PAPER=de_BE.UTF-8          LC_NAME=de_BE.UTF-8
 [9] LC_ADDRESS=de_BE.UTF-8        LC_TELEPHONE=de_BE.UTF-8
[11] LC_MEASUREMENT=de_BE.UTF-8    LC_IDENTIFICATION=de_BE.UTF-8

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets
[8] methods   base

other attached packages:
 [1] httr_1.4.1                  stringr_1.4.0
 [3] readr_1.3.1                 xlsx_0.6.3
 [5] VariantAnnotation_1.32.0    Rsamtools_2.2.3
 [7] Biostrings_2.54.0           XVector_0.26.0
 [9] SummarizedExperiment_1.16.1 DelayedArray_0.12.2
[11] BiocParallel_1.20.1         matrixStats_0.55.0
[13] Biobase_2.46.0              GenomicRanges_1.38.0
[15] GenomeInfoDb_1.22.0         IRanges_2.20.2
[17] S4Vectors_0.24.3            BiocGenerics_0.32.0
[19] optparse_1.6.4

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.2               lattice_0.20-38          prettyunits_1.0.2
 [4] xlsxjars_0.6.1           assertthat_0.2.1         zeallot_0.1.0
 [7] digest_0.6.15            BiocFileCache_1.10.2     R6_2.2.2
[10] backports_1.1.4          RSQLite_2.1.1            pillar_1.4.2
[13] zlibbioc_1.32.0          rlang_0.4.0              GenomicFeatures_1.38.2
[16] progress_1.2.2           curl_4.2                 blob_1.1.1
[19] Matrix_1.2-17            RCurl_1.98-1.1           bit_1.1-14
[22] biomaRt_2.42.0           compiler_3.6.1           rtracklayer_1.46.0
[25] pkgconfig_2.0.3          askpass_1.1              openssl_1.4.1
[28] tidyselect_0.2.5         tibble_2.1.3             GenomeInfoDbData_1.2.2
[31] XML_3.98-1.11            crayon_1.3.4             dplyr_0.8.3
[34] dbplyr_1.4.2             GenomicAlignments_1.22.1 bitops_1.0-6
[37] rappdirs_0.3.1           grid_3.6.1               DBI_1.0.0
[40] magrittr_1.5             stringi_1.4.3            getopt_1.20.3
[43] vctrs_0.2.0              tools_3.6.1              bit64_0.9-7
[46] BSgenome_1.54.0          glue_1.3.1               purrr_0.3.2
[49] hms_0.5.1                AnnotationDbi_1.48.0     memoise_1.1.0
[52] rJava_0.9-11

Thats a coders issue, I'd contact the authors


text:
Topmed File Requirements
I am trying to run my VCF's through Topmed Imputation (https://imputation.biodatacatalyst.nhlbi.nih.gov/#!pages/home)
I am trying to make sure that I have the correct requirements for my file to be able to run through the software. These requirements are:
Create a separate vcf.gz file for each chromosome.
Variations must be sorted by genomic position.
GRCh37 or GRCh38 coordinates are required.
I have managed to split my VCF into each chromosome but I am unsure how to check which coordinates are being used and if my file is sorted by genomic position. Are there unix commands I can use to check, or do I just look in the file, Ive looked and not noticed anything.

There are tools that allows you to sort your VCF file by genomic position like bcftools sort or vcf-sort. But if you don't want to install these software, you could always do with bash tools for one vcf:
# get the header of your VCF and write to file 'sorted.vcf'
# if your file has the extension .gz, use zgrep instead
grep "^#" your.vcf > sorted.vcf

# get the rest of the lines not starting with '#', sort them first by chromosome, then position and add them to file 'sorted.vcf'
grep -v "^#" your.vcf | sort -V -k1,1 -k2,2n >> sorted.vcf

But since you have one file for each chromosome and assuming that your files look like sample_chrZ.vcf.gz, where Z is the chromosome number, you will have to do a loop:
mkdir sorted_vcf
for chr in *.vcf.gz; do
   filename="${chr%.*}"
   zgrep "^#" "$chr" > sorted_vcf/"$filename"
   zgrep -v "^#" "$chr" | sort -V -k1,1 -k2,2n  >> sorted_vcf/"$filename"
   bgzip sorted_vcf/"$filename"
done

And if you want to know whether your VCF has been aligned against GRCh38 or GRCh37, check the header in the VCF file with a line starting with ##reference=. A command-line like grep "##reference=" your.vcf should give you the result. However, if you don't have the reference mentioned in your VCF file, this is going to be harder to guess. This post should help you figure it out.

Hi, are you familiar with TopMed?

Hello, no I am not, sorry. I answered your questions according to your need.


text:
connect same ids with values in id rows (result of effectiveELD)
I can't cope with reworking the result. If id is the same I want to stick those items.
with open(x,'r') as rf:
    l = fr.readlines()
    previous_line_id = ''
    corrected_lines = []

    for i in l:
       i_split = i.split('\t')
       id_ = i_split[0]
       if id_ == previous_line
           xxxxxxxxxx
       else:
           corrected_lines .append(i)

result:
    WP_XXXXXXXX2.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5
    WP_192XXXXXX.1  PF13489 Methyltransferase domain    4
    XXXXXXX45.1 PF02310 B12 binding domain  12  <-----------  DOUBLE IDS
    XXXXXXX45.1 PF04055 Radical SAM superfamily 5 <-----------
    WP_XXXXXXX29.1  PF13489 Methyltransferase domain    4
    WP_0XXXXXXXX.1  PF13489 Methyltransferase domain    4
    WP_XXXXXXX90.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5  <--------------   DOUBLE ID
    WP_XXXXXXX90.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  7  <--------------
    WP_0XXXXXX62.1  PF13489 Methyltransferase domain    4

result I need:
WP_XXXXXXXX2.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5
WP_192XXXXXX.1  PF13489 Methyltransferase domain    4
XXXXXXX45.1 PF02310 B12 binding domain  12; PF04055 Radical SAM superfamily 5 <----------- STICKED SAME ID,OTHER VALUES
WP_XXXXXXX29.1  PF13489 Methyltransferase domain    4
WP_0XXXXXXXX.1  PF13489 Methyltransferase domain    4
WP_XXXXXXX90.1  PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5; PF06325  Ribosomal protein L11 methyltransferase (PrmA)  7  <--------- STICKED STICKED SAME ID,OTHER VALUES
WP_0XXXXXX62.1  PF13489 Methyltransferase domain    4

I am sorry, but I don't understand what you want exactly. Could you please elaborate on what you are trying to achieve ?

@user324810 I commented text

@user324810 I want to remove the double id and merge their values together

Thank you. Give me a moment

Assuming that the file test.txt is tab-separated values with these fields:
WP_XXXXXXXX2.1  PF06325     Ribosomal protein L11 methyltransferase (PrmA)      5
WP_192XXXXXX.1  PF13489     Methyltransferase domain                            4
XXXXXXX45.1     PF02310     B12 binding domain                                  12
XXXXXXX45.1     PF04055     Radical SAM superfamily                             5
WP_XXXXXXX29.1  PF13489     Methyltransferase domain                            4
WP_0XXXXXXXX.1  PF13489     Methyltransferase domain                            4
WP_XXXXXXX90.1  PF06325     Ribosomal protein L11 methyltransferase (PrmA)      5
WP_XXXXXXX90.1  PF06325     Ribosomal protein L11 methyltransferase (PrmA)      7
WP_0XXXXXX62.1  PF13489     Methyltransferase domain                            4

You can proceed first by reading your file:
with open("test.txt",'r') as rf:
    l = rf.readlines() # all the lines are stored in variable 'l'

You should create a dictionary with your ID as keys and the remaining line as values:
# create empty dictionary
dict_id = {}

for i in l:                                                               # for each line
    key_id = i.split("\t")[0]                                             # get the id
    if key_id not in dict_id.keys():                                      # if id as key is not found in dictionary
        value = "\t".join(i.split("\t")[1:]).rstrip()                     # parse the value
        dict_id[key_id] = [value]                                         # add the id as key to dictionary with its value
    else:                                                                 # else if id as key is already in dictionary
        dict_id[key_id].append("\t".join(i.split("\t")[1:]).rstrip())     # append to id as key

And now, you can print the result:
for key, value in dict_id.items():         # for each key, values in dictionary
    print(key + ": " + "; ".join(value))   # print the key and their values separated by '; '

And you should obtain:
WP_XXXXXXXX2.1: PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5
WP_192XXXXXX.1: PF13489 Methyltransferase domain        4
XXXXXXX45.1: PF02310    B12 binding domain      12; PF04055     Radical SAM superfamily 5
WP_XXXXXXX29.1: PF13489 Methyltransferase domain        4
WP_0XXXXXXXX.1: PF13489 Methyltransferase domain        4
WP_XXXXXXX90.1: PF06325 Ribosomal protein L11 methyltransferase (PrmA)  5; PF06325      Ribosomal protein L11 methyltransferase (PrmA)  7
WP_0XXXXXX62.1: PF13489 Methyltransferase domain        4

The final script:
with open("test.txt",'r') as rf:
    l = rf.readlines()

dict_id = {}
for i in l:
    key_id = i.split("\t")[0]
    if key_id not in dict_id.keys():
        value = "\t".join(i.split("\t")[1:]).rstrip()
        dict_id[key_id] = [value]
    else:
        dict_id[key_id].append("\t".join(i.split("\t")[1:]).rstrip())

for key, value in dict_id.items():
    print(key + ": " + "; ".join(value))

realy thanks :) all best


text:
Is there a publication database or search engine offering geneID or UniProtID correlation?
Do we have a publication database or search engine that offers geneID or UniProtID correlation?
For example, I search "IPS cells" -> I want to know what genes or proteins are under active research in IPS cells.

I think this could be close to what you're looking for: http://glad4u.zhang-lab.org/index.php
Paper: GLAD4U: deriving and prioritizing gene lists from PubMed literature
It works with gene names though, so I am not sure it's a 100% fit.

Take a look at GeneWalk, https://churchman.med.harvard.edu/genewalk. It's a very interesting solution to this rather difficult problem. It is well documented and packaged. But It's open source software -- not a service. You'll have to run it on your own system.


text:
Extracting a column from a vcf file
I have a zipped vcf file of dbSNP hg38 version
No space left on my device to unzip that
I want to extract a column from that
I have tried this
cut -f 3 00-All.vcf.gz | sed 's/[\t]/,/g' > cols.csv

But that needs my vcf to be unzipped
Do you know any solution please?

Either the zcat command or gunzip -c will decompress the file in a streaming fashion to the standard output without writing the decompressed data to disk. Pipe the output to cut, sed, etc.
Try
zcat 00-All.vcf.gz | cut -f 3 | sed 's/[\t]/,/g' > cols.csv

or
gunzip -c 00-All.vcf.gz | cut -f 3 | sed 's/[\t]/,/g' > cols.csv

But if you just want a single 'column' (I presume you mean field), then you are better off just doing bcftools query.

On some systems I have issues with zcat but gunzip -c will also accomplish the same result.

@Exhausted Don't use unix utilities for this - it needs a level of expertise that your post history suggests you do not have. Use bcftools instead, that will address a lot of edge cases that you can't foresee.

tar can do select zip file search and extraction, but you're also data munging.

Thank you all, but why you down voting me?

I have No idea, I didn’t downvote.