text:
How critical is the "Non-unique read group integer IDs found" warning in running ipdSummary?
I am trying to get the ipdRatio and modified bases information from a PacBio bam file which is sorted and indexed. The command I'm running is:
ipdSummary sorted_chr21_0_100000.bam --reference $REF_path/$REF_FILE_NAME.fasta --referenceWindow 20:0-100000 --identify m6A,m4C,m5C_TET --methylFraction --outfile sorted_chr21_0_100000;

I get the following warning message:
RuntimeWarning: Non-unique read group integer IDs found - some features may be restricted warnings.warn("Non-unique read group integer IDs found - some features may be restricted", RuntimeWarning)

My question is which types of features might be restricted if my read groups are not unique?
I am not sure how critical it could be.


text:
Trouble Running VEP in Conda Environment Due to Perl DBD::mysql Error
Hello StackExchange Community,
I am currently facing issues with running the Variant Effect Predictor (VEP) and vcf2maf in a Conda environment due to a Perl-related error.
Installation Steps: I installed VEP and vcf2maf using Conda with the following commands, and everything installed without any problems:
conda create -n vcf2maf -y python=3.8
conda activate vcf2maf
conda install -y  -c bioconda vcf2maf
conda install  -y  -c bioconda ensembl-vep
conda install  -y   -c bioconda samtools
vcf2maf.pl  --help
vep --help
samtools -v
perl -e '{print join"\n",@INC}

Cache Installation: I then installed the cache for VEP using these commands:
export VEP_PATH=$HOME/vep
export VEP_DATA=$HOME/.vep
mkdir -p $VEP_PATH $VEP_DATA
cd $VEP_PATH
vep_install -a cf -s homo_sapiens -y GRCh38 -c  $VEP_PATH --CONVERT
#  downloading ftp://ftp.ensembl.org/pub/release-88/variation/VEP/homo_sapiens_vep_88_GRCh38.tar.gz
cd $HOME/.vep nohup wget
ftp://ftp.ensembl.org/pub/release-101/variation/indexed_vep_cache/homo_sapiens_vep_101_GRCh38.tar.gz &
tar zxvf homo_sapiens_vep_101_GRCh38.tar.gz

However, when I tried to test run VEP with the following command:
vep --help
# By default the VEP uses $HOME/.vep/
vep  -i test_varscan.vcf  -o test.vcf  \
--cache --dir_cache $HOME/vep --force_overwrite  --assembly GRCh38 --vcf

I encountered the following error:
perl: relocation error: /public2/qiuting/.conda/envs/vcf2maf/lib/perl5/site_perl/5.22.0/x86_64-linux-thread-multi/auto/DBD/mysql/mysql.so:symbol mysql_options4, version libmysqlclient_18 not defined in file libmysqlclient.so.18 with link time reference

To resolve this, I tried:
cpanm DBD::mysql
conda install -c conda-forge perl-dbd-mysql

Unfortunately, these attempts did not resolve the issue, and the same error occurred.
I am not very familiar with Perl and its modules. Could anyone please help me solve this issue? Any suggestions or insights would be greatly appreciated!
Thank you in advance!

More details about the configurations/ version should be needed.  Ting, pls use the following question as a reference:  https://stackoverflow.com/questions/60987782/facing-issue-while-configuring-mysql-with-apache-airflow-in-hadoop

source: Stack Exchange

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

sort: votes


text:
ipdSummary command stops working if --identify m5C_TET is included in command, and does not give any m5C info either
I have installed the latest SMRT Tools (ipdSummary version 3.0). My goal is to extract ipdRatio and modified bases. In the ipdSummary --help command :
--identify IDENTIFY   Specific modifications to identify (comma-separated list). Current options are m6A, m4C, m5C_TET. Using --control overrides this option. (default: m6A,m4C)

Accordingly, I am running this command:
ipdSummary sorted_chr21_000.bam --reference $REF_path/$REF_FILE_NAME.fasta --referenceWindow 20:1-200000 --identify m6A,m4C,m5C_TET -j 20 --methylFraction --outfile $save_path/$chrom/$name

My bam file is sorted and relatively small (<200 MB).
But this command never finishes and does not throw any errors. In different stages of working it just stops. I checked different thread numbers and smaller files, but still, the command doesn't work.
However, once I changed 'm5C_TET' to 'm5C' it worked perfectly fine and very fast on the same bam file, but in the output, I could not see any information on 'm5C' modifications. I checked the new command with many other bam files and still the same.
My question is the following:
How can I get at least 'm5C' modifications? Should I change my SMRT Tools version?

I asked the same question on the ipdSummary GitHub page and I got the following response:
"I think that was an experimental option that was never fully working or validated, and it no longer works. We probably need to update the help text. Use the separate tool jasmine if you want to do m5C calling."
jasmine also is a tool in the SMRT toolkit that analyzes the kinetic signatures of cytosine bases in CpG motifs to identify the presence of 5mC.

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


text:
Printing the results of for loop into a txt file
I have a for loop which prints available vcf files in the path for a list of patient IDs (using find dx data function) into my zsh terminal:
for i in $(cat patient_id.txt); do
 dx find data --property patient_id=$i --path "vcf/"
done

How I can have a file contain the vcf files that the loop shows in terminal? I mean by using this command in the terminal, I am seeing the name and path of my vcf files for a given participant list and working directory; Instead of seeing the file name in the terminal, I want to have that file name in a txt file.
Thanks

Please [edit] your question and explain what happens if you simply add > file after the done. Is output redirected to file? Is it still shown in your terminal? What is the output? Are you sure files are found and you aren't just getting error messages?

This command worked just by adding  > file after done
for i in $(cat patient_id.txt); do
dx find data --property eid=$i --path 'vcf/'
done > vcf_list_file.txt

Things that write to standard output on the terminal can be wrapped up in brackets and redirected to a file. Assuming this is what you want to do, I added brackets around the existing code, then used > to redirect the output to a file.
(for i in $(cat patient_id.txt); do
 dx find data --property patient_id=$i --path 'vcf/'
done) > vcf_files.txt

Note that the brackets go around the entire code, including the for bit.

@zizogolu the brackets go around the entire code, including the for bit

@Zizogolu can you please type out what you put into the terminal, and the error messages that were produced?

There is no need to make this run in a subshell (that's what the parentheses do). All you need is the redirect in the end. for i in ...; do ...; done > file works and redirects all stdout of the for loop. Compare for i in a b; do echo "$i"; done > file and the same thing with parentheses, it makes no difference other than running the whole thing in a subshell.

Thanks a million @terdon , your comment solved my issue

That is... odd, @Zizogolu. If it works with for ...; do ... done > file then it would also work with Gringer's (for ...; do ... done) > file, the parentheses are just irrelevant, they wouldn't stop the redirection from working. In any case, can you please post what worked for you as an answer so the question can be marked as answered?

I suspect the issue was in the path command-line argument, which seems to have changed since I copied the code. I've updated that now.


text:
Featurecount output .txt file from bam file
After running this command:
featureCounts -p -s 2 -a $genome -o $dir"/"$specie"/count_table/"$value".txt" $output_loc -T 8

on the BAM output of HISAT (output_loc), I got this output:

Now, does the last column of this file show the gene count value?
And if I merge the same column of different samples, how should I go through the normalization step?
Thanks a lot for your help!

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

The first 6 column in standard featureCounts output represent what is in the column names. All columns after than (starting at 7) represent the counts for the sample(s). If you use a single bam file as input then it's one column, if you use many bams as input then it is one column per bam. If the counts are gene level, exon transcript or anything else depends on how you run the tool. Typically, with default settings, it is usually the raw gene level counts.
Generally, you can sum up technical replicates (=lane/sequencing replicates to increase depth) prior to running normalization.

Thank you so much I appreciate your help!
My other question is that since I didn't know that I could give multiple bam files to the feature count I instead ran the feature count on each bam file and then merged column 7 of the output text files and made my final count table in which rows are genes and columns are samples and cells are the raw gene counts. Now, as I know I should do normalization at this step, would TPM be a good method? if it's so off the table can you please let me know the sources I should use to learn how to proceed?
Thanks a lot!

You simple list all bam files to be processed.
featureCounts (options...) one.bam two.bam / or just *.bam
Read manuals of edgeR or DESeq2, use these normalizations. TPM is poor. Many posts on this at biostars.org why that is.

Thanks, ATpoint!


text:
How do I get the gene annotations as a text file from a genbank file?
Also posted on biostars
I have a genbank file and I need the following details about the gene annotations in a text file
# ID    name    start   stop    frame   product description function    cogcats status_region   status_function type    region_number   function_number

How do I go about extracting them? Please let me know if more information is needed from me! Thank you

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

Please [edit] your question and add i) an example input genbank file and ii) the exact output you would want from that example. Also tell us what operating system you are using so we know what tools you have available. Will this file describe a transcript? A gene? A chromosome? Most of the fields you mention aren't standard, so we really need the example to understand.

The fact that Genbank files are also text files, I believe you meant that you prefer the same data in the "table format". In that sense why not fetching your annotations of interest in GFF/GTF format? And if that is not available for your genome of interest, there are tools that convert Genbank to GFF. Seeing the column names in your example, I think you can go for the GFF3 format and "expand" the last column (attritubes) to extract the semicolon-separated info within that column.

Crossposted here: https://www.biostars.org/p/9591298/


text:
I'm trying to run aTRAM tool for assembly, but i'm stuck on this error:
Code:
$ python3 atram.py --max-processes=12 \
  --blast-db=/dgxb_home/se23plsc009/Project/results/aTRAM_db/T1.sqlite.db \
  --query-split=/dgxb_home/se23plsc009/Project/data/ref_orthologs/2395_query_file \
  --output-prefix=/dgxb_home/se23plsc009/Project/results/Temporary/Output.fasta-unitigs.fa \
  --assembler=abyss \
  --path=/dgxb_home/se23plsc009/softwares/aTRAM/lib/assemblers/velvet.py \
  --iterations=3 \
  --max-target-seqs=3000 \
  --log-file=/dgxb_home/se23plsc009/Project/results/aTRAM_assembly1/log_file \
  --temp-dir=/dgxb_home/se23plsc009/Project/results/aTRAM_assembly1/Temporary/ \
  --keep-temp-dir \
  --protein

Output:

2024-03-29 08:21:29.100271 INFO : aTRAM blast DB = "/dgxb_home/se23plsc009/Thrips_Project/results/aTRAM_db/T1", query = "2395_query_file_RPRC015459_PA_2395.fasta", iteration 1
2024-03-29 08:21:29.101075 INFO : Blasting query against shards: iteration 1
2024-03-29 08:23:00.952156 INFO : All 93 blast results completed
2024-03-29 08:23:00.954933 INFO : 72 blast hits in iteration 1
2024-03-29 08:23:00.954965 INFO : Writing assembler input files: iteration 1
2024-03-29 08:23:01.258918 INFO : Assembling shards with abyss: iteration 1
2024-03-29 08:23:01.406976 ERROR: Exception: [Errno 2] No such file or directory: '/dgxb_home/se23plsc009/Thrips_Project/results/aTRAM_assembly/Temporary/atram_nowsbxk1/T1_2395_query_file_RPRC015459_PA_2395.fasta_01_clkdfnle/output.fasta-unitigs.fa'
2024-03-29 08:23:01.409907 INFO : Writing 0 filtered contigs after iteration 1
2024-03-29 08:23:01.412235 INFO : 0 total contigs after iteration 1

I have tried running the job with different assembler, gave the path to "Output.fasta-unitigs.fa" as output folder but nothing is working.
I've attached below the github link of the tool for reference: https://github.com/juliema/aTRAM/blob/master/doc/atram_reference.md
I have checked to make sure that each output directory exists.
Here is the github issue that I think is related to my problem:
https://github.com/juliema/aTRAM/issues/308


text:
Dante Labs VCF: UPDATE: SOLVED
Has anyone tried to analyse the VCF files provided by Dante Labs? (Human WGS, Illumina)
I'm looking into .raw.vcf and .filtered.vcf.
For reference sequence, all I can see is 2 contigs:
GL000225.1 and GL000192.1.
All variants belong to one of these.
There's no conventional CHROM column.
I don't know how those coordinates map to chromosomes in GRCh37 or GRCh38 assemblies.

GL000192 corresponds to chr1_gl000192_random
GL000225 corresponds to chrUn_gl000225

Moreover, there's no VCF header with metainfo, all I can see is 1 top row
:11,27:0.711:38:4,11:7,16:25:64,0,24:2.8811e+01,1.4126e-02,2.7139e+01:0.00,34.77,37.77:6,5,12,15:7,4,18,9.

What's the use of those VCFs?
What variants do they contain?
What version of VCF it is and where can I get the standard metainfo?
Is there a map of those coordinates to any core assembly
chromosomes?
Am I right they contain not all genomic VCFs but only a small part
of the genome? (Size of GL000225.1 contig is 211,173 bp and size of
GL000192.1 is 547,496 bp, that obviously doesn't sum up to the whole
genome).

Many thanks for your replies!
UPD: I found a mistake, sorry. The file was incomplete, only tail was uploaded to Jupyterhub.
In local Jupyter, the file loaded completely, it actually has the  header and CHROM column.


text:
Trying to plot the results of GSEA in python
I'm having some major chaos with the output table generated and then passing that to my plot function. New to this (1st time), could somebody review this piece of code and suggest some corrections. I saved the output of a data frame as a txt file and it does not look good at all. The numbers and columns are there, but they follow some sort of strange wa, like a staircase.
# Create DataFrame from significant results
df_go = pd.DataFrame(
    [(
        r.GO,
        r.goterm.name,
        r.goterm.namespace,
        r.p_uncorrected,
        len(r.study_items),  # Number of study genes
        GO_items.count(r.GO),
        ', '.join([rev_mapper[y] for y in r.study_items])  # Concatenate study genes into a single string
    ) for r in goea_results_sig],
    columns=['GO', 'term', 'p', 'n_genes', 'n_study', 'n_go', 'study_genes']
)

print(df_go)
# Save DataFrame to a text file (CSV format)
df_go.to_csv('df_go_results.txt', sep='\t', index=False)

I can provide more code, because when I tried to commit and push from IDE onto github I got the following message.
Successfully created project 'GSEAplots' on GitHub, but initial commit failed: inflate: data stream error (incorrect data check) corrupt loose object 'a87848b96903648398f723df040fd0be540a025d' unable to read a87848b96903648398f723df040fd0be540a025d
No idea why or what ? Git does not seem to interact with PyCharm and/ or Github really well ?
The plots code :
import matplotlib.pyplot as plt
import numpy as np

plt.figure(figsize=(6,4),dpi=300)
ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

plt.rcParams["font.size"] = 18
plt.rcParams["font.family"] = "Helvetica"

plt.scatter(df_p['Z'],-np.log10(df_p['P(sgRNAs)']),s=18,color='grey',alpha=0.3,linewidths=0.3,edgecolors='grey')

label_ls1 = list(df_up['Gene'])

df_show1 = df_p.loc[label_ls1,:]
lfc_ls1 = list(df_show1['Z']); p_ls1 = list(-np.log10(df_show1['P(sgRNAs)']))

plt.scatter(lfc_ls1,p_ls1,s=20,color='red',linewidths=0.5,edgecolors='grey')

for i, txt in enumerate(label_ls1[-4:]):
    plt.annotate(txt, (lfc_ls1[-4:][i]+0.03, p_ls1[-4:][i]-0.08),fontsize=8)

label_ls2 = list(df_down['Gene'])

df_show2 = df_p.loc[label_ls2,:]
lfc_ls2 = list(df_show2['Z']); p_ls2 = list(-np.log10(df_show2['P(sgRNAs)']))

plt.scatter(lfc_ls2,p_ls2,s=20,color='blue',linewidths=0.5,edgecolors='grey')

for i, txt in enumerate(label_ls2[1:5]):
    plt.annotate(txt, (lfc_ls2[1:5][i]+0.03, p_ls2[1:5][i]-0.08),fontsize=8)

plt.axhline(y=1.3, linestyle='--', color='lightgray')
plt.axvline(x=-2, linestyle='--', color='lightgray')
plt.axvline(x=2, linestyle='--', color='lightgray')

plt.xlabel('Differential score (PRMT5/DMSO)')
plt.ylabel('-log10(p-value)')

#plt.savefig('DifferentialEssential_volcanoplot.png',dpi=300,facecolor='white')
plt.show()

Data frame that I am trying to plot has the following structure :
             total_samples  occurrence_total  occurrence_CRSsNP  \

ENSG00000100592             15                15                  8
ENSG00000101336             15                14                  8
ENSG00000134285             15                15                  8
             occurrence_CTL  sparsity  rlog_av_CTL  rlog_av_CRSsNP  \

ENSG00000100592               7  0.099193    11.709476       11.413235
ENSG00000101336               6  0.322330     4.084709        6.179892
ENSG00000134285               7  0.098323     9.169283        9.595925
            rlog_Regulation  vst_av_CTL  vst_av_CRSsNP vst_Regulation  \

ENSG00000100592             CTL   11.790292      11.367062            CTL
ENSG00000101336          CRSsNP    4.909465       7.051789         CRSsNP
ENSG00000134285          CRSsNP    9.142557       9.734598         CRSsNP
              raw_av_CTL  raw_av_CRSsNP Raw_Regulation structural_zero  \

ENSG00000100592  3973.285714       2428.500            CTL              No
ENSG00000101336    10.285714        119.750         CRSsNP              No
ENSG00000134285   602.857143        761.625         CRSsNP              No
             DESeq2Norm_av_CTL  DESeq2Norm_av_CRSsNP DESeq_Regulation  \

ENSG00000100592        3526.244455           2614.486202              CTL
ENSG00000101336           9.293777            131.796045           CRSsNP
ENSG00000134285         535.929110            828.762618           CRSsNP
                baseMean  log2FoldChange     lfcSE      stat  \

ENSG00000100592  3039.973387       -0.431423  0.076683 -5.626030
ENSG00000101336    74.628320        3.828859  0.605862  6.319689
ENSG00000134285   692.106981        0.628415  0.105894  5.934376
             pvalue_DESEq2  padj_DESEq2        pvalue      padj  \

ENSG00000100592   1.844041e-08     0.000120  3.694635e-06  0.038127
ENSG00000101336   2.620902e-10     0.000005  2.009572e-07  0.006221
ENSG00000134285   2.949646e-09     0.000029  1.052822e-06  0.016297
             errorBarUp  errorBarLow Final_Regulation external_gene_name

ENSG00000100592   -0.354739    -0.508106              CTL              DAAM1
ENSG00000101336    4.434721     3.222997           CRSsNP                HCK
ENSG00000134285    0.734309     0.522521           CRSsNP             FKBP11

Hi @aleksk31, welcome. It’s quite hard to answer your question as there in fact multiple questions there, and not enough context. Try and think how to break down your problem into smaller parts. If you are having GitHub or IDE issues open a new question on the main stack exchange. I suggest editing your question down to just your first immediate question, which is about saving a pandas dataframe. When you edit, tell us where the variable goea_results_sig comes from and why the resulting csv isn’t what you expect (we don’t know what you mean by staircase)


text:
Generating CIGAR strings for MSA with reference
I have a list containing the results of a MSA:
msa = [
    '-G-ATG-CGGATACCG-', # reference sequence
    'CGTATCTCG--TACCGA',
    '-GTATCTCGT-GACCG-',
    'CG-AT--CGG-GACCGA',
    '---AT--CG--GACCGA',
    '-GTAT---GC-GACCG-',
    'CG-AT--CGG-GACCGA',
    'CGTAT--CG--GA CG-',
]

I would like to generate the CIGAR strings of those alignments (against the reference sequence); the expected result would be something like:
cigars = [
    '13M',
    '1I1M1I3M1I2M2D5M1I',
    '1M1I3M1I3M1D5M',
    '1I3M1D3M1D5M1I',
    '1D2M1D2M2D5M1I',
    '1M1I2M2D2M1D5M',
    '1I3M1D3M1D5M1I',
    '1I1M1I2M1D2M2D5M',
]

I don't know what would be a good approach for this. I managed to get the "aligned pairs" of each alignment but I don't think that'll help for determining the cigar strings:
# mapping = position in alignment => position in reference
mapping = {y: x for x,y in enumerate(y for y,nuc in enumerate(msa[0]) if nuc != '-')}

for aln in msa:
    aligned_pairs = [(y, mapping[z]) for y,z in enumerate(z for z,nuc in enumerate(aln) if nuc != '-') if z in mapping]
    print(aligned_pairs)

[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (7, 7), (8, 8), (9, 9), (10, 10), (11, 11), (12, 12)]
[(1, 0), (3, 1), (4, 2), (5, 3), (7, 4), (8, 5), (9, 8), (10, 9), (11, 10), (12, 11), (13, 12)]
[(0, 0), (2, 1), (3, 2), (4, 3), (6, 4), (7, 5), (8, 6), (9, 8), (10, 9), (11, 10), (12, 11), (13, 12)]
[(1, 0), (2, 1), (3, 2), (4, 4), (5, 5), (6, 6), (7, 8), (8, 9), (9, 10), (10, 11), (11, 12)]
[(0, 1), (1, 2), (2, 4), (3, 5), (4, 8), (5, 9), (6, 10), (7, 11), (8, 12)]
[(0, 0), (2, 1), (3, 2), (4, 5), (5, 6), (6, 8), (7, 9), (8, 10), (9, 11), (10, 12)]
[(1, 0), (2, 1), (3, 2), (4, 4), (5, 5), (6, 6), (7, 8), (8, 9), (9, 10), (10, 11), (11, 12)]
[(1, 0), (3, 1), (4, 2), (5, 4), (6, 5), (7, 8), (8, 9), (9, 10), (10, 11), (11, 12)]

How can I generate the CIGAR strings efficiently?

If you loop through the characters of two aligned sequences in parallel then you can easily determine the operations that compose the CIGAR. Then you just have to "collapse" the generated list of operations:
import itertools

ref = msa[0]

for aln in msa:
    cig = []

    for i,nuc in enumerate(aln):
        if nuc != '-':
            if ref[i] != '-':
                cig.append('M')
            else:
                cig.append('I')
        elif ref[i] != '-':
                cig.append('D')

    cig = [(len(list(it)), op) for op,it in itertools.groupby(cig)]

    print( ''.join(str(n) + op for n,op in cig) )

13M
1I1M1I3M1I2M2D5M1I
1M1I3M1I3M1D5M
1I3M1D3M1D5M1I
1D2M1D2M2D5M1I
1M1I2M2D2M1D5M
1I3M1D3M1D5M1I
1I1M1I2M1D2M2D5M


text:
How does BioJava retrieve bond information from a PDB file?
How can I compute energy from a protein PDB file?
I asked the above question.
Then I thought, Why don't I check how BioJava obtains bond information from a protein PDB file?
Surprisingly, BioJava's tutorial doesn't have any discussion on Bonds.
I checked BioJava's GitHub repo, but I failed to find any such information.
How does BioJava retrieve bond information from a PDB file?


text:
Ramachandran Plot for specific amino acids in a pdb file
Does anyone know of a Ramachandran plot server (URL please) that allows you to display, for example, just the leucine residues in a protein structure rather than displaying all the amino acids in the pdb file? Thanks!
Paul

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

https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/midas/ramachandran.html  : Usage: ramachandran [ assign ]  atom-spec
The ramachandran command shows the distribution of peptide φ,ψ angles (a Ramachandran plot) for each specified model that contains protein. It is the command-line implementation of the Ramachandran plot function in the Model Panel. A separate plot will be generated for each protein-containing model with any atom in atom-spec, and the plot for a model will contain all of its peptide φ,ψ angles regardless of which atoms were specified.


text:
Parameter optimisation of minimap2
I'm leveraging minimap2 to map select genes from short-read microbial fastq metagenomics zip files.
minimap2 --heap-sort=yes --splice-flank=no metagenome.fastq.gz fastaDB.fas

What I noticed is the low tolerance of mismatch when I need to maximise mismatch because microbes have loads of genetic diversity. Thus the collection of references sequences are never really a reference, in fact the "microbe reference genome" is a misnomer, because in situ diversity of microbes is huge.
Question How do you relax the mapping stringency criteria for minimap2? If you are familiar with Blast parameters this would help alot here. I know blast but I've never parameterised it.
The options are this are:
Alignment options

-A Matching score [2]

-B INT    Mismatching penalty [4]

-O INT1[,INT2] Gap open penalty [4,24].

-s INT    Minimal peak DP alignment score to output [40].

--score-N INT
Score of a mismatch involving ambiguous bases 1

Anything else that is relevant.
Note, -f is a famous option, not sure it helps here.
I'm going to code up a hyperparameterisation unit-test, but understanding the parameters first would really, really help.
Alternative approaches welcome as well.

Plugging minimap2. I code high performance (hadoop-style, back-tracking, graph, traverse tall buildings in a stride  etc ...). The coder that wrote minimap has knowledge beyond, i.e. very fast run times directly from gzip file. I assume they are indexing then chunking the raw zip file then unzipping each 'chunk' JIT (just in time).

Thanks @MaximilianPress great answer. I'll fill in the backstory.
I'm trapping total genetic diversity across Genbank and the metagenomes. My protocol to do both:

Blasting out Genbank and then
piping the output across metagenomes via in this case minimap2.

Hence speed is key. So I ain't replacing Blast with minimap2 but leveraging minimap2 via the output of blast.
However on reflection of @MaximilianPress's discussion the answer is to generate a gold standard of diversity to assess alternative approaches.

First of all, use short read option
The short read preset (-x sr) for minimap2 is heavily recommended for doing such mapping:
minimap2 -x sr
Heng Li on minimap2 for related genomes
Heng Li (minimap2 author) has some suggestions on a somewhat related github issue (though it is now 4 years old):

Minimap2 is not really intended for this purpose. For short reads he recommends stampy, but it is clearly quite old (python2.7) and does not seem to have much recent developer history.
Failing that, he suggests an option set such as follows:

If you want to give minimap2 a try anyway, you can use the following setting:  minimap2 -axsr -A2 -B4 -O4,24 -E2,1
If minimap2 is not sensitive enough, reduce -k and -w at the cost of mapping speed. Reducing -w should have relatively lower impact on speed, but reducing -k should be more sensitive.

I would suggest engaging with Heng, he is quite responsive in general and may have more suggestions.
Other aligners for metagenomics
I would suggest that you look for more possibilities for alignment, or possibly not using alignment in the first place. However, this is very sensitive to your ultimate application, which I'm not sure that I understand from your post.
I was pleased to see that others have reported good results using minimap2 for metagenomics, which I was skeptical about, having been out of the field for a bit. Usually in metagenomics there is an expected tradeoff between sensitivity and specificity for the reasons that you state (this paper is old but it uses BLAST, so it is a sort of best case for sensitivity on unrelated genomes at the level of divergence that they test).
There are of course new shotgun metagenomic aligners coming out all the time, here is one recent example that I know nothing about but which may be interesting as a comparison.
If you are really looking for a rich parameter set that allows low-stringency mappings, I would suggest going back to BLAST, though it will be much slower. minimap2 and related algorithms are to some extent dependent on the ability to find short seed matches, so there will necessarily be less than perfect sensitivity in your application.

Many thanks. Cool answer. I've addressed the issues you've below the final "-------"

Just to follow-up on the performance of minimap2: the solution @MaximillianPress described works. Thousands of hits are obtained for comparatively obscure genes (some will be duplicates). Parameterisation gives a 2-log10 increase (astonishing).
The issue is that its:

relatively slow - precisely as described by @MP;
very, very sensitive to parameterisation, just a single integer changes the output by a couple of thousand hits.

I can see why there was concern about using it - such notable sensitive is an issue - and alternatively its clear why other investigators have had great success with it - hit the "sweet-spot" and its a bonanza. For my purposes it's fine - it's not a quantitation of complete genetic diversity, as long as I capture a reasonable proportion of diversity that is fine. Clearly, they'll be times the "sweet-spot" is missed but thats still okay (famous last words :-) ).


text:
BLAST+ not working on my Command Line
As you know, ncbi has several databases, including nt, nr, etc.
What I am going to do is to blast a sequence of approximately 30 base pairs to the nr database and extract all similar sequences.
For this purpose, with the help of YouTube, I downloaded and installed Command Line Blast on my Windows PC.
Then, from the database section I downloaded all the nr databases, which were 83 zip files and 200 GB in size.  Then, I unzipped it and stored it to a folder on my computer.
Because I saw that the size of the first two zip files that came out of the zip increased from 25 GB to 50 GB, I only unzipped two files to see if my command line works or not.
So my database became a folder of almost 50 GB.
Then I ran the simplest command in the cmd of my computer.
blastn -query D:\query.fa -out D:\results.txt -db D:\nrDataBase

It's been almost two weeks now that I'm stuck day and night with this simple line of code and it doesn't work.
It says once that the query file is not accessible. Once it says that no suitable files were found in the database, etc.
I tested every solution that came to my mind, I used the text format for the query, but it didn't work
I downloaded one or two files from the nt database, I said maybe there is a problem with the nr database, but it still didn't work.
I made a small fasta file, introduced it as a database, but it still didn't accept.
I tried also on the Macbook, it still didn't work, either the query got an error or the database.
The blast has been well installed, When using -h it shows the full info of the tool. It doesn’t work means that I encounter different errors like: The query is not acceptable Or The alias or index file is found in the database I encounter these errors while e. g., the query file is a normal fasta file. I mean, the errors don’t seem to be logical, it looks like something else is going wrong.
I'm really upset, I told you to send me a message to see if you can help me.
Thank you very much and sorry for the long message.

"it doesn't work" is not enough information. What do you see on the command prompt when you just run blastn -h?

@M__ the question was migrated; the OP does not have a profile on bioinfo SE yet.

Thanks @RamRS of course I should have spotted.

You state that you are not using the entire database, you're using just a subset of the resource files. That could be causing the error. Just because you're using a subset of the resource files does not mean you're picking a functional subset of the database.

Please provide more context around what you want to do. It is not necessary to download the entirety of the nr database in order to extract only the sequences that match a specific subsequence. There's probably a better way to achieve your goal without massive downloads.

try running makeblastdb on your downloaded database files, then re-run your blast command

Are you working in the correct directory? Meaning, are you in the directory of the script you are trying to run? What did the error message look like? Please run again and provide that detail so it gives us more to work with!:)


text:
Nanopore Flongle vs MinION
How do the Oxford Nanopore Flongle flow cells perform compared to the standard MinION flow cells, when it comes to reliability and error rates? I found this comparison article from last year, from where the below graph is taken, but they used the old chemistry. How is the situation today with the new R10.4.1 nanopores and Kit 14 chemistry? Does anyone have first-hand experience, or a good reference?
My use case would mainly be metabarcoding, i.e. amplicon sequencing. If anyone has experience of Flongle flow cells for specifically metabarcoding, I would be happy to hear your thoughts.

For us, R10.4.1 Flongle flow cells are running at reduced yield due to lower numbers of working pores (something like 10-20 pores instead of 30-60 pores for R9.4.1 flow cells). This doesn't really affect the stuff we're using Flongles for, because I make sure we've got plenty of headroom built in when designing experiments.
Here are the yield assumptions I use for flow cells, which are basically the lower ranges that we've seen when using them:

Flongle - 200 Mb
MinION - 5 Gb
PromethION - 50 Gb

For Flongle, this has been fine for us to do sequencing and assembly of 12 plasmids on one flow cell using the rapid barcoding kit, and the improved accuracy from the Kit14 + R10.4.1 combination makes reliable assembly easier (especially where homopolymers are involved). I expect it will work similarly well for amplicon sequencing.
Error rates are pretty similar across platforms, because it's the pore that matters, not the package. Here's a comparison plot of mapped accuracy vs predicted accuracy (from the basecaller) that I did recently on a TrackIt 1kb Plus DNA ladder, with bead cleanup using LFB, run on a PromethION flow cell with the V14 ligation kit (sequenced on a P2 Solo), then called using the new bacterial model. The reads represented in this plot have been filtered to only keep (at most) the longest 2000 matches for each target. The run produced 1,128,671 reads (1.95 Gb) over 1.5 hours:

Here's a comparison plot of mapped accuracy vs predicted accuracy (from the basecaller) for the same DNA ladder sample run on a Flongle flow cell (sequenced on a MinION Mk1b), with bead cleanup using SFB, also taking at most 2000 longest sequences from each mapping. The run produced 465,144 reads (191 Mb) over 24 hours:

The predicted accuracy is taken from the basecaller (i.e. fastq quality string), and is calculated by summing up the quality values after converting to a linear representation, then dividing by the number of bases and converting back to a log representation:
sub getMeanQual {
  my ($qual) = @_;
  if(length($qual) == 0){
    return(0);
  }
  my $qBase = 33;
  my @qspl = split("", $qual);
  my $qualTotal = 0;
  my $qualCount = 0;
  foreach my $qual (@qspl) {
    $qualTotal += 10 ** (-(ord($qual) - $qBase) / 10);
    $qualCount++;
  }
  return(int(-10 * log($qualTotal / $qualCount) / log(10)));
}

According to the trained models that LAST generated, there's not much different between the called reads between Flongle and PromethION flow cells. Here's the PromethION base mismatch matrix (q45 SNPs):
# substitution percent identity: 99.9971
#last -Q 1
#last -t4.38507
#last -a 23
#last -A 19
#last -b 2
#last -B 4
#last -S 1
# score matrix (query letters = columns, reference letters = rows):
       A      C      G      T
A      6    -45    -35    -46
C    -45      6    -45    -44
G    -33    -45      6    -46
T    -46    -45    -46      6

And here's the Flongle base mismatch matrix (q41 SNPs):
# substitution percent identity: 99.9922
#last -Q 1
#last -t4.36483
#last -a 20
#last -A 20
#last -b 4
#last -B 5
#last -S 1
# score matrix (query letters = columns, reference letters = rows):
       A      C      G      T
A      6    -45    -36    -46
C    -45      6    -44    -45
G    -32    -45      6    -46
T    -46    -45    -46      6

So I guess you could say that the Flongle has slightly lower accuracy, but I wouldn't call that difference significant for any practical purposes.

Thank you for the update, I have marked your answer as a solution!


text:
Compare files with similar names using placeholders/wildcards in Snakemake
In Snakemake I want to use two input files for comparison. The one file has a sample name (e.g. 526846), with an experimental id (e.g. EF_S13' and 'MF_S5'). The file I want to compare it to has the same sample name only. I want to compare each sample with experiment ID (first files) to the file with corresponding sample name (second file). These files are in different directories. I thought I would be able to use two expand functions with two lists, as shown in the code below.
Note: these are example lists and the lists in reality is muchlarger.
Relevant code:

SAMPLE_SET = ['12345', '56789', '54321']
ID_SET = ['EF_S13', 'MF_S5', 'EF_S11', 'MF_S3', 'EF_S10', 'MF_S2']

rule all:
    input:
        array = expand("results/{sample}_{samp_id}_array/{sample}_{samp_id}.summary.csv", zip, sample = SAMPLE_SET, samp_id = ID_SET)

rule happy_rtg_array:
    """ Normalizes and compares vcf files with array samples"""
    input:
        query_vcf = expand("WESP_run/results/variant_calling/{sample}_{samp_id}.vcf", zip, sample = SAMPLE_SET, samp_id = ID_SET),
        truth_vcf = expand(os.path.join(config["truth_dir_vcf_gsa"],"{sample}.vcf.gz"), sample = SAMPLE_SET),
    output:
        "WESPValidation/results/{sample}_{samp_id}_array/{sample}_{samp_id}.summary.csv"
[more code]

However, when I look at the output log file with --dry-run, I get as input for the job a huge list with all the inputs in one go, rather then per file/comparison, and one output file (see relevant output), and that times the amount of the (right) combinations of the file name and id (counted from the output).
Also the wildcards are not what I thought they should be. I might be confused with the expand function, as I thought how I used them now should "lock in" the wildcards and become placeholders filled in by the lists.
Relevant output:
[Tue Oct  3 14:31:20 2023]
rule happy_rtg_array:
    input: 12345_EF_S13.vcf, 12345_MF_S5.vcf, 56789_EF_S11.vcf, 56789_MF_S3.vcf, 54321_EF_S10.vcf, 54321_MF_S2.vcf, 12345.vcf, 56789.vcf, 54321.vcf
    output: results/12345_EF_S13_array/12345_EF_S13.summary.csv
    log: happy_rtg_array/12345_EF_S13.log
    jobid: 6
    reason: Missing output files: WESPValidation/results/12345_EF_S13_array/12345_EF_S13.summary.csv
    wildcards: sample=12345_EF, samp_id=S13
    threads: 10
    resources: mem_mb=10000, disk_mb=7855, tmpdir=/tmp, threads=1, runtime_min=120

I doubt this is how it should be, and wonder what I am doing wrong. Is there a better/correct way to do this.
Thank you in advance

After some trial and error, I fixed it with the expand function only in the rule all and removing the expand from the input in the rule happy_rtg_array (thanks to @Pallie), and by changing the output file paths:
SAMPLE_SET = ['12345', '56789', '54321']
ID_SET = ['EF_S13', 'MF_S5', 'EF_S11', 'MF_S3', 'EF_S10', 'MF_S2']

rule all:
    input:
        array = expand("results/{sample}/{samp_id}/summary.csv", zip, sample = SAMPLE_SET, samp_id = ID_SET)

rule happy_rtg_array:
    """ Normalizes and compares vcf files with array samples"""
    input:
        query_vcf = "_".join(["WESP_run/results/variant_calling/{sample}","{samp_id}.vcf"]),
        truth_vcf = os.path.join(config["truth_dir_vcf_gsa"],"{sample}.vcf.gz"),
    output:
        "WESPValidation/results/{sample}/{samp_id}/summary.csv"
[more code]

Basically, I replaced all the underscores with slashes. For some reason these underscores are messing up the wildcards when trying to connect them with it. Maybe someone more knowledgeable knows why this is.
[Tue Oct  5 12:24:10 2023]
rule happy_rtg_array:
    input: 12345_EF_S13.vcf
    output: results/12345/EF_S13/summary.csv
    log: happy_rtg_array/12345_EF_S13.log
    jobid: 6
    reason: Missing output files: WESPValidation/results/12345/EF_S13/summary.csv
    wildcards: sample=12345, samp_id=EF_S13
    threads: 10
    resources: mem_mb=10000, disk_mb=7855, tmpdir=/tmp, threads=1, runtime_min=120

For the reason why, check the part about ambiguous wildcards: https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#wildcards. I missed the amibguity of your wildcards when I wrote my answer :)

Dank je @Pallie :-)

Using expand in your rule's input:  causes snakemake to aggregate all vcfs into one list and leads to your undesired behaviour. By simply removing the expand calls from the input: you should get your desired result.
https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#the-expand-function

This is only partially the answer, I tried this before indeed, but still get the following wildcards "wildcards: sample=12345_EF, samp_id=S13" and they should be: "wildcards: sample=12345, samp_id=EF_S13", to get the right file name for the input "truth_vcf"


text:
ValueError: count not convert string to float
I'm trying to use a Python script to plot with a rolling window using pandas and seaborn. This code worked for the longest time but now it's giving me an error that I don't know how to fix.
Here is the code I'm using:
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

heterozygosity_df = pd.read_csv("file.tsv", sep="\t", header=None).set_axis(['chr', 'pos', 'het'], axis=1, copy=False)

heterozygosity_df['het'] = heterozygosity_df['het'].transform(lambda x: -1*(x-1))

windows = [25, 50, 100, 200]
steps = [5]

for window in windows:
    for step in steps:
        heterozygosity_ratio = {}
        for chr_heterozygosity_df in heterozygosity_df.groupby('chr'):
            heterozygosity_ratio[chr_heterozygosity_df[0]] = chr_heterozygosity_df[1].set_index('pos').rolling(window=window, step=step, center=True).mean()
...

There is more after this but that's just to configure the resulting plot with matplotlib. The error I get when I use the above code is this:
Traceback (most recent call last):
  File "/opt/homebrew/lib/python3.11/site-packages/pandas/core/window/rolling.py", line 368, in _prep_values
    values = ensure_float64(values)
             ^^^^^^^^^^^^^^^^^^^^^^
  File "pandas/_libs/algos_common_helper.pxi", line 42, in pandas._libs.algos.ensure_float64
ValueError: could not convert string to float: 'chrI'

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/opt/homebrew/lib/python3.11/site-packages/pandas/core/window/rolling.py", line 485, in _apply_blockwise
    arr = self._prep_values(arr)
          ^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/homebrew/lib/python3.11/site-packages/pandas/core/window/rolling.py", line 370, in _prep_values
    raise TypeError(f"cannot handle this type -> {values.dtype}") from err
TypeError: cannot handle this type -> object

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "plot_heterozygosity_window.py", line 29, in <module>
    heterozygosity_ratio[chr_heterozygosity_df[0]] = chr_heterozygosity_df[1].set_index('pos').rolling(window=window, step=step, center=True).mean()
                                                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/homebrew/lib/python3.11/site-packages/pandas/core/window/rolling.py", line 2111, in mean
    return super().mean(
           ^^^^^^^^^^^^^
  File "/opt/homebrew/lib/python3.11/site-packages/pandas/core/window/rolling.py", line 1507, in mean
    return self._apply(window_func, name="mean", numeric_only=numeric_only)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/homebrew/lib/python3.11/site-packages/pandas/core/window/rolling.py", line 615, in _apply
    return self._apply_blockwise(homogeneous_func, name, numeric_only)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/homebrew/lib/python3.11/site-packages/pandas/core/window/rolling.py", line 487, in _apply_blockwise
    raise DataError(
pandas.errors.DataError: Cannot aggregate non-numeric type: object

I have tried converting the 'chr' column directly from a string to a float:
heterozygosity_df['chr'] = heterozygosity_df['chr'].astype(float)
I have tried converting the entire dataframe to a float:
heterozygosity_df = heterozygosity_df.astype(float, errors = 'raise')
But everything keeps giving me the ValueError: could not convert string to float: 'chrI'
Does someone know an easy way to fix this? Like I said, I have been using this code for months, and it's always worked - have they deprecated something in the pandas package to make it so you can't read in columns with strings and integers?
EDIT:
For reference the dataframe is structured like this:
chrI  123  1
chrI  234  0
chrI  237  1
chrII 456  1
chrII 567  1
chrII 678  1
...
chrXVI 123 1

Where the first column designates the chromosome, the second column is the position, and the third column is the heterozygosity indicator.

The following would be useful,
print (heterozygosity_df.columns)

This looks like static typing issue within pandas. What I suspect is heterozygosity_df.['chrI'] is a column in the dataframe. What I think has happened is there's a mix of strings and floats within this column. pandas has set this as a "string" but you are wanting to perform numerical operations. Thus the solution is simply
print(heterozygosity_df.dtypes) # this should state 'chrI' is a "category" or "string"
heterozygosity_df['chrI'] = heterozygosity_df['chrI'].astype(float)
print(heterozygosity_df.dtypes)

If you have multiple changes the syntax is
heterozygosity_df = heterozygosity_df.astype({'chrI':'float', 'egColumn':'category'})

I suspect there will be other errors, e.g. the header is the first row of the data column. This is because again for pandas to automatically assign a column to a "string" means there must be a string value within the column.

From the comments. I see whats happening. The easy solution is ...
replacement = {
    "chrI": 1,
    "chrII": 2,
    "chrIII": 3,
    "chrIV": 4,
    ....
}

heterozygosity_df['chr'] = heterozygosity_df['chr'].str.replace(replacement, regex=True)

From the comments ... good it works ... this is what I would have personally done ..
replacement = {
    "chrI": 1,
    "chrII": 2,
    "chrIII": 3,
    "chrIV": 4, # continue for all chromosomes
}

heterozygosity_df = pd.read_csv("file.tsv", sep="\t", header=None).set_axis(['chr', 'pos', 'het'], axis=1, copy=False)
heterozygosity_df['chr'] = heterozygosity_df['chr'].str.replace(replacement, regex=True).astype('int')

Normally str should be in place because when 'chr' is imported its an object. However, if it works thats the only thing that counts.

when I do print(heterozygosity_df.columns) it outputs: Index(['chr', 'pos', 'het'], dtype='object') and when I do print(heterozygosity_df.dtypes) it outputs: chr object    pos  int64    het  int64    dtype: object

Since I have chromosomes going from chrI to chrXVI I grouped the chromosomes by "chr" instead of individual "chrI" so when I did heterozygosity_df['chr'] = heterozygosity_df['chr'].astype(float) it just gave me the same ValueError: could not convert string to float: 'chrI'

the annoying part is that I'm not using the 'chr' column to do the actual numerical operations - I just use it to group the data in the beginning (see: for chr_heterozygosity_df in heterozygosity_df.groupby('chr'):) the actual math goes into the position and heterozygosity columns

@rimo Addressed below the "-----"

Traceback (most recent call last):
File "fixing_error.py", line 90, in
heterozygosity_df['chr'] = heterozygosity_df['chr'].str.replace(replacement, regex=True)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/opt/homebrew/lib/python3.11/site-packages/pandas/core/strings/accessor.py", line 129, in wrapper
return func(self, *args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: StringMethods.replace() missing 1 required positional argument: 'repl'

I get this error message when I implement the code you wrote above @M__

wait I got it to kind of work by taking out the '.str' and leaving the command like this: heterozygosity_df['chr'] = heterozygosity_df['chr'].replace(replacement, regex=True)

however now it's only printing chromosomes 1, 5, 10, and 17 (chrmt) for whatever reason and not the rest of them???

Ah! I finally got it to work! I used your command but I had to change it up a little bit to work with my data. The final command I used is this:
replacement = {
    "chrI": 1,
    "chrII": 2,
    "chrIII": 3,
    "chrIV": 4,
    "chrV": 5,
    "chrVI": 6,
    "chrVII": 7,
    "chrVIII": 8,
    "chrIX": 9,
    "chrX": 10,
    "chrXI": 11,
    "chrXII": 12,
    "chrXIII": 13,
    "chrXIV": 14,
    "chrXV": 15,
    "chrXVI": 16,
    "chrmt": 17
}

heterozygosity_df['chr'] = heterozygosity_df['chr'].replace(replacement, regex=False)

And with that I was able to generate plots showing all of the chromosomes! (I have to figure out how to change the axis labels from 1,2,3,4...10 back to the chromosome but that's a future problem). Thank you so much for all your help!!

regex=False you might be right.


text:
Standardizing Variant IDs best practise
I'm currently pondering on the best approach to handle and standardize variant IDs within our department to mitigate the challenges associated with rsids and the potential ambiguity of just using chromosomal positions.
The GWAS glue (https://mrcieu.github.io/gwasglue2/articles/Strategy.html) packagae, state oa method where a variant ID is crafted by integrating the chromosomal position (chr), effect allele (ea), and non-effect allele (nea). It assumes the format: chr:pos_ea_nea. For alleles that exceed 10 characters, like certain indels, they apply the murmur32 hashing algorithm from the R digest package.
Given our objectives:
We likely prefer not to have excessively lengthy IDs.
Yet, using hashed IDs for certain long variants might be perplexing for some who aren't informed about the approach.
Hashing every ID could yield consistency, but perhaps might not be as intuitive for regular use.
I'm keen on collecting insights:
Thoughts on the described method from Bristol?
What might be the maximum length for an indel we could encounter?
Is there a broader effort underway to address this concern?
Would adopting the Bristol approach set a potential standard, or should we explore or develop alternative methods? For example should the genome build be built into the id? or species?
I appreciate the community's feedback as we strive for clarity and efficiency in our variant naming conventions.
For clarity:
I want to store data on variants. For example GWAS summary statistics.
I want to be able to match variants from differnt studies so I want to have an identifier for the variants.
I could use rs numbers
I could use chr:pos
I could use chr:pos:ref:alt
I could use chr:pos:effect:noneffect
I could use chr:pos:population:effect:noneffect
or studyid:chr:pos:populatoin:effect:noneffect
or etc
There are pros and cons for each of these but I would like to here what people on here think these are.
I am also intrested in is there work out there to do this already. e.g. to update the idea of rsids to consider multiple alleles and effects in populations?

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

I have tried to clarify

I would not hash. Variant IDs should be static across builds, similarly to how rsIDs are handled between b37 and b38. I do not think hashing is worth it because there's not that much downside for a long ID

You need this to work on multiple species, right? And across different reference genome builds? Are you only interested 8n species with a sequenced genome or do you also need to handle variants defined with respect to a specific transcript or protein only, in the absence of a reference genome?

Also, are you open to paid tools? The company I work for has an API that will, among other things, provide unique ids for arbitrary (human) variants. If you are open to non-free tools check out https://landing.varsome.com/varsome-api.

Thanks ill check that out.

Thoughts on the described method from Bristol?

I don't immediately see the benefit to this over rsIDs.

What might be the maximum length for an indel we could encounter?

Potentially very long. However, why does it matter? This will be rare.

Is there a broader effort underway to address this concern?

dbSNP comes to mind.

Would adopting the Bristol approach set a potential standard, or
should we explore or develop alternative methods? For example should
the genome build be built into the id? or species?

It depends on the type of analysis you're planning on doing. Are you doing cross-species analyses? And if so, what types? Do you expect to be comparing genomes with different genome builds?

For example GWAS summary statistics. I want to be able to match
variants from different studies so I want to have an identifier for
the variants.

Both chr:pos as well as rsID work for this, but just make sure that the genome build matches if you go with chr:pos.

I could use rs numbers I could use chr:pos I could use chr:pos:ref:alt
I could use chr:pos:effect:noneffect I could use
chr:pos:population:effect:noneffect or
studyid:chr:pos:populatoin:effect:noneffect or etc

If you're sure your input data has rs numbers and always will, that might be easiest. However, you may find that sometimes it doesn't have rs numbers (e.g. some Affymetrix data), and it's surpassingly annoying to add them. I'd choose chr:pos, personally, and then consider the effect/non-effect/ref/alt/ on a per-analysis basis, because I would want to still identify the same position if there was a different ref/alt for example.
Edit:
Based on the comment (OP wants to match effect direction), I would suggest chr:pos:effect_allele, but of course this means there is an effect allele (i.e. the position has been studied).

The vast majority of variants (remember that variants are effectively  infinite since you can have any length of indel at any position) won't have rsIDs. Novel variants won't, by definition. So if the OP needs to work with arbitrary variants, rsIDs won't help. If not, rsIDs are indeed the best approach.

@terdon Yes. I think chr:pos would be ideal/most simple, and OP would just have to be careful regarding builds

Thanks for the response and your thoughts.
I guess the thing is in the context of merging GWAS results for coloc and mr etc, chr:pos does not seem sufficeint. I can make sure everything in on the same build, but how do I know if an effect of a snp on say bmi in a japanese study can be matched to a  study which finds an effect at the same position in a european study. I need to take into account ref/alt to effect non effect , direction..

@user27815 I see. I would try chr:pos:effect_allele in that case. Would that work?

gnomAD uses chr-pos-ref-alt, not as an official ID but as a unique way to access web pages detailing each variant. That might be worth looking into. Example: https://gnomad.broadinstitute.org/variant/1-55051215-G-GA?dataset=gnomad_r3

chr-pos-ref-alt is a good way to match variants, but may become problematic (without further tinkering) if the same variant has opposite effect sizes in different studies

yes, which is why I think some times chr:pos:effect:noneffect is often used, then that clears that up. But you still have the issue of population if the effect is differnt in differnt pops..


text:
Error with eggNOG-mapper jobs
I recently tried using eggNOG-mapper for pan-genome and core genome analysis of files from OrthoFinder, specifically the 'Orthogroup_Sequences' which contains a huge list of protein files in fasta format. As instructed, the folder was compressed into a tar.gz zip file and uploaded. If the upload was successful, an email is sent to your mail ID where you can check the details of your job and decide whether to START the job or not. Every time I START the job, it fails within a few minutes. Initially I presumed that it was due to the file size and hence reduced the number of sequences submitted, but in vain. Can anyone help me?
eggNOG-mapper : http://eggnog-mapper.embl.de/

Which version of eggNOG-mapper is being used here?

Hi @M__ I think it's v2. That's what the citation says at the bottom of the page. it's probably v2.1.12. Do you think the command line version would be better?

There is a known bug between eggNOG-mapper2 and eggNOG6. This is described here:
https://github.com/eggnogdb/eggnog-mapper/issues/428
In your case, I'm not sure on what strategy you are deploying after OrthoFinder, but there's a chance it's not solvable at present and the best place for information is the eggnog GitHub.

Thank you @M__ . I have posted my issue on the github page for failing eggNOG-mapper jobs. I have also tried to perform the same function using COGclassifer after concatenating all the .fa sequences in the Orthogroup_sequences file. Hopefully it yields reliable results. If you have any other suggestions please let me know.

@K_081 Thanks it's a positive step forward. eggNOG is very, very good. I'm not be looking for a solution outside this. I'll keep an eye out for your GitHub post and I will be looking around for a way forward and post back here once done (timeline several weeks).

Hi @M__, I just got a reply on the github page from one of the developers I guess, and he has found an unwanted line in my file.  It is the file created by OrthoFinder itself so I don't know how it appeared there. I'll have to check all the files before uploading them to eggNOG-mapper.

I'd like to answer my own question as I have found the solution for this and thought it might help anyone else trying to use eggNOG-mapper.
The prerequisite for uploading the input file with the required collection of .faa files, is that, the file extension has to end with .gz. The site mentions that too. I had made the mistake of uploading a .tar.gz which is what my system produced as a compressed file of my folder of sequences. Although the upload was successful, the runs failed as the file had some unwanted lines owing to the .tar.gz format.
I randomly uploaded a concatenated .faa file which worked well and got the results too! The trick is that, the input file has to be strictly a .gz and NOT .tar.gz or it could also be a .faa file with all the sequences.

Thanks for this that is very helpful.


text:
How to interpret a GenomeScope plot?
I am working on a new sponge genome and I have produced the genomescope plot from the 21-mer kmer frequencies. I have hard time interpreting the plot. Can someone please help me?
Thank you.

Hi @MudithEkanayake IMO the question has been answered. Please ask via comments if you're not clear about the answer. Do remember to upvote/accept when you're happy with the answer.

Based on the numbers shown above the plot, and the appearance of the curve, this looks like reads from a diploid genome approximately 113Mb in length that was sequenced at around 174X coverage.


text:
Compare substitution matrices
I am looking for methods to compare substitution matrices to determine which one is better for a given set of data / type of protein.
Do you know such methods / tools?
e.g.  BLOSUM62 vs PHAT for 50k - 100k+ sequences

The technique is phylogenetics is called likelihood ratio and for a regular amino acid matrix is 2Δln which is assessed on Chi2 tables where Δln is difference in likelihood between model 1 and model 2. The degrees of freedom hmmmm... is the number of parameters that are constrained, this involves the number of taxa minus the number of parameter differences between the models. This makes performing 50 000 taxa ... tricky because the degrees of freedom will be outside the Chi2 table.
iqtree will perform a tree and provide a likelihood for 50+K taxa. There are two limitations,

No PHAT model
degrees of freedom (as mentioned)

Overcoming df is simply to say whatever is the lowest likelihood is the best model and simply state Δln without bothering with significance testing. If the difference between models is large you'll get away with that (needs to be big) and for 50+K taxa it will be need to be huge.
Getting PHAT (the model below the ---- might be PHAT) into iqtree is certainly doable because of QMaker. However, how this is done I don't know, because I never done this calculation with iqtree, you'll need to read the manual. For traditional packages such as PAML it's easy ... you make a tree then specific the custom amino acid model and it spits out the likelihood. For the number of taxa involved, PAML is not possible, it's impossibly big for PAML, even 1000 is a lot of PAML.
iqtree does have LG, WAG and JTT amino acid models and will almost certainly have all the BLOSUM matrices and you just specify those at the command line. Again a comparison with PHAT needs to specify a custom matrix and I don't know how to do that in iqtree. The only thing I know is it must be doable.
iqtree is the most efficient of the maximum likelihood phylogenetic methods by a very long way. Hence for dealing with 50+K is why I'm recommending it here.
Note The theory about whats going on and why the whole thing needs phylogenetics is very complicated. Basically its an entire lecture on phylogenetic theory.

This I believe is PHAT (needs checking)
https://github.com/biopython/biopython/pull/283
Note the code is copied and pasted from the above link, therefore it should be strictly regarded as a quote for copyright purposes.
#!/usr/bin/env python
matrix = {
('W', 'F') :  0.0, ('L', 'R') : -6.0, ('I', 'I') :  5.0, ('Q', 'Q') :  9.0,
('W', 'N') : -5.0, ('V', 'I') :  3.0, ('H', 'T') : -4.0, ('H', 'P') : -6.0,
('W', 'V') : -4.0, ('Q', 'E') :  1.0, ('W', 'R') : -7.0, ('Q', 'A') : -3.0,
('H', 'H') : 11.0, ('H', 'D') : -1.0, ('L', 'N') : -3.0, ('Y', 'M') : -2.0,
('Y', 'I') : -3.0, ('Y', 'E') : -2.0, ('E', 'S') : -3.0, ('Y', 'A') : -3.0,
('Y', 'Y') : 11.0, ('T', 'C') : -1.0, ('E', 'C') : -7.0, ('Y', 'Q') :  0.0,
('E', 'G') : -3.0, ('V', 'A') :  1.0, ('C', 'C') :  7.0, ('M', 'R') : -6.0,
('P', 'T') : -4.0, ('V', 'E') : -5.0, ('P', 'P') : 13.0, ('I', 'T') : -1.0,
('K', 'S') : -5.0, ('R', 'G') : -5.0, ('I', 'P') : -4.0, ('R', 'C') : -8.0,
('A', 'T') :  0.0, ('K', 'K') :  5.0, ('A', 'P') :  3.0, ('V', 'M') :  1.0,
('I', 'D') : -5.0, ('K', 'C') :-10.0, ('K', 'G') : -5.0, ('R', 'S') : -6.0,
('F', 'Q') : -2.0, ('F', 'A') : -1.0, ('V', 'V') :  4.0, ('M', 'N') : -2.0,
('F', 'E') : -5.0, ('D', 'N') :  2.0, ('F', 'I') :  0.0, ('F', 'M') :  0.0,
('M', 'S') : -2.0, ('S', 'S') :  6.0, ('L', 'Q') : -3.0, ('W', 'E') : -7.0,
('W', 'A') : -4.0, ('W', 'M') : -4.0, ('H', 'S') : -2.0, ('W', 'I') : -4.0,
('S', 'C') :  1.0, ('L', 'A') : -1.0, ('L', 'E') : -5.0, ('W', 'Q') :  1.0,
('H', 'G') : -4.0, ('Q', 'N') :  2.0, ('H', 'C') : -7.0, ('L', 'M') :  2.0,
('W', 'Y') :  1.0, ('Y', 'N') :  2.0, ('E', 'P') : -5.0, ('Y', 'F') :  4.0,
('E', 'T') : -5.0, ('A', 'A') :  5.0, ('I', 'N') : -3.0, ('G', 'A') :  1.0,
('Y', 'V') : -3.0, ('E', 'D') :  6.0, ('W', 'H') : -3.0, ('Y', 'R') : -6.0,
('M', 'Q') : -1.0, ('P', 'S') : -3.0, ('R', 'H') : -4.0, ('A', 'C') :  1.0,
('R', 'D') : -7.0, ('K', 'P') : -4.0, ('L', 'D') : -5.0, ('K', 'T') : -6.0,
('V', 'N') : -3.0, ('M', 'A') : -1.0, ('K', 'H') : -5.0, ('V', 'R') : -7.0,
('P', 'C') : -8.0, ('M', 'E') : -5.0, ('A', 'S') :  2.0, ('T', 'T') :  3.0,
('R', 'T') : -6.0, ('I', 'G') : -2.0, ('R', 'P') : -7.0, ('K', 'D') : -5.0,
('I', 'C') : -3.0, ('F', 'R') : -7.0, ('F', 'V') : -1.0, ('L', 'C') : -2.0,
('F', 'F') :  6.0, ('D', 'A') : -5.0, ('F', 'N') : -1.0, ('W', 'D') : -7.0,
('L', 'P') : -5.0, ('Q', 'S') : -1.0, ('N', 'C') : -2.0, ('N', 'G') : -1.0,
('H', 'N') :  4.0, ('W', 'T') : -7.0, ('Q', 'G') : -2.0, ('W', 'P') : -6.0,
('Q', 'C') : -5.0, ('N', 'S') :  1.0, ('L', 'H') : -4.0, ('L', 'L') :  4.0,
('G', 'T') : -1.0, ('M', 'M') :  6.0, ('G', 'P') : -3.0, ('Y', 'K') : -4.0,
('Y', 'G') : -3.0, ('Y', 'C') : -1.0, ('E', 'A') : -5.0, ('E', 'E') : 12.0,
('Y', 'S') : -2.0, ('M', 'P') : -5.0, ('V', 'C') : -2.0, ('M', 'T') :  0.0,
('V', 'G') : -2.0, ('R', 'E') : -6.0, ('V', 'K') : -8.0, ('K', 'Q') : -1.0,
('R', 'A') : -6.0, ('I', 'R') : -6.0, ('N', 'A') : -2.0, ('V', 'S') : -2.0,
('M', 'D') : -5.0, ('M', 'H') : -4.0, ('K', 'A') : -7.0, ('R', 'Q') : -2.0,
('K', 'E') : -4.0, ('F', 'S') : -2.0, ('I', 'K') : -7.0, ('D', 'P') : -5.0,
('D', 'T') : -5.0, ('I', 'M') :  3.0, ('F', 'C') :  0.0, ('W', 'L') : -3.0,
('F', 'G') : -2.0, ('F', 'K') : -7.0, ('F', 'T') : -2.0, ('D', 'D') : 12.0,
('Q', 'T') : -3.0, ('W', 'G') : -5.0, ('Q', 'P') : -3.0, ('W', 'C') : -4.0,
('W', 'K') : -8.0, ('H', 'Q') :  2.0, ('Q', 'D') :  0.0, ('W', 'W') : 11.0,
('V', 'L') :  1.0, ('L', 'G') : -2.0, ('W', 'S') : -5.0, ('L', 'K') : -7.0,
('N', 'P') : -4.0, ('H', 'E') : -1.0, ('N', 'T') : -1.0, ('H', 'A') : -3.0,
('Y', 'L') : -2.0, ('Y', 'H') :  3.0, ('G', 'S') :  1.0, ('Y', 'D') : -4.0,
('V', 'Q') : -3.0, ('L', 'T') : -1.0, ('G', 'G') :  9.0, ('G', 'C') : -2.0,
('E', 'N') :  0.0, ('Y', 'T') : -3.0, ('Y', 'P') : -5.0, ('R', 'N') : -3.0,
('V', 'D') : -5.0, ('K', 'R') : -1.0, ('V', 'H') : -5.0, ('I', 'Q') : -3.0,
('V', 'P') : -4.0, ('M', 'C') : -2.0, ('K', 'N') : -2.0, ('V', 'T') :  0.0,
('M', 'G') : -1.0, ('T', 'S') :  1.0, ('I', 'E') : -5.0, ('M', 'K') : -6.0,
('I', 'A') :  0.0, ('N', 'N') : 11.0, ('R', 'R') :  9.0, ('F', 'P') : -5.0,
('L', 'I') :  2.0, ('I', 'S') : -2.0, ('D', 'S') : -4.0, ('L', 'S') : -2.0,
('I', 'H') : -5.0, ('F', 'D') : -5.0, ('D', 'C') : -7.0, ('F', 'H') : -2.0,
('D', 'G') : -2.0, ('F', 'L') :  1.0
}

chars = [aa[0] for aa in matrix] + [aa[1] for aa in matrix]
chars = sorted(list(set(chars)))

output = " ,%s\n" % ",".join(chars)

for i, aa1 in enumerate(chars):
    output += aa1
    for j, aa2 in enumerate(chars):
        if (aa1, aa2) in matrix:
            output += "," + str(matrix[(aa1, aa2)])
        else:
            output += "," + str(matrix[(aa2, aa1)])
    output += "\n"
print(output)


text:
TCGA pathology slides magnification
Does anybody know if there is a way to find out what are the magnification levels of the TCGA pathology slides? Are they all 40x? I can't seem to find a label of that anywhere

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

Can you provide a link to any images? It's likely that the resolution is available in the file metadata.

Thanks. Here is an example of one of them. I can't seem to find the resolution in the metadata

https://portal.gdc.cancer.gov/files/f23db5f3-6b2e-471e-9b18-e477140d9a5e


text:
How to fix "No gene can be mapped, Expected input gene ID, return NULL..." error in clusterProfiler?
I am trying EnrichKEGG on a few gene samples of genes to check it works before inputting the full genome:
This is my data input as test.csv
id,fc
AJAP_14870,-0.04
AJAP_14875,0.32
AJAP_14880,-0.06
AJAP_14885,0.15
AJAP_14890,-0.52
AJAP_14895,-0.39
AJAP_14900,0.16

geneList = test[,2]
names(geneList) = as.character(test[,1])
geneList = sort(geneList, decreasing = TRUE)

mygene <- names(geneList)[abs(geneList) > 2]

kk <- enrichKEGG(gene         = mygene,
                 organism     = 'aja',
                 pvalueCutoff = 0.05)

There are no spaces in the gene names, I tested via
print(head(test$ID), quote=TRUE)

But I get the error
--> No gene can be mapped....
--> Expected input gene ID:
--> return NULL...

I also tried installing through other means
getOption("clusterProfiler.download.method")
#[1] "libcurl"
install.packages("R.utils")
library(R.utils)
R.utils::setOption("clusterProfiler.download.method","auto")

but that gives me the same error when I run the code. Any ideas how to fix this?
Thanks

Your code produces no error on my machine, what is your clusterProfiler version ?

@Basti Thanks so much for getting back to me!!

packageVersion("clusterProfiler")
[1] ‘4.2.2’

Maybe you could try BiocManager::install("clusterProfiler", force = TRUE)as mine is ‘4.7.1’ and restart your R session

@Basti Thanks for this! I tried that but it still just says version 4.2.2??

Could you share your sessionInfo() in your post ? Maybe your R version or Bioconductor version is incompatible with latest clusterProfiler version.

@Basti you're absolutely right - I've just had to update rtools, biocmanager and R! It's working now! Thanks so much for your help!! The relief haha

No worries, I hope you will get good results!

@Basti, would you mind converting your comments to an answer? It would be more visible.

This issue may arise with earlier versions of clusterProfiler.
A solution is to re-install clusterProfiler : BiocManager::install("clusterProfiler", force = TRUE)
Pay attention that updating your R or Bioconductor (and restart your R session) version may be needed as noted in the comments

enrichKEGG() expects entrez gene IDs so you will need to convert your gene identifiers to entrez identifiers first.
EDIT: Please see the answer by @Basti.

I've used this code previously and it has only worked using the locus_tags as thats what is input into KEGG for aja as far as I'm aware? Really confused why its not working now?


text:
Where to download a file with major and minor alleles at every position?
I want a list of all variants,  i.e. sites which are known to vary between human to human. For example, it should ideally cover all sites in here, but without samples.
I don't want a giant reference panel with many samples, nor a .fa which has no major/minor allele information.
I just want a file with the following information at each locus:

Location (chromosome position and chromosome)
Bonus (not necessary): rsID
Major allele (i.e., A, C, T, or G) at that position
Minor allele at that position

I'm not picky about how major/minor alleles are defined. Usually they are defined as more common in the general population vs. less common (loosely). Conforming to b37 or b38 would be ideal.
I'm looking for one such file that doesn't have samples in it, but does contain all possible positions (not all VCF files have this)
Where can I download something like this? I like working with VCFs but any format is good.

Sounds like you're looking for the dbSNp files. You haven't specified a species, but you can download the data for human from the dbSNP FTP site. Specifically, you want:
wget https://ftp.ncbi.nlm.nih.gov/snp/organisms/human_9606/VCF/00-All.vcf.gz

This file has lines like:
1   10848   rs1271898990    C   G   .   .   RS=1271898990;RSPOS=10848;dbSNPBuildID=151;SSR=0;SAO=0;VP=0x050000020005000002000100;GENEINFO=DDX11L1:100287102;WGT=1;VC=SNV;R5;ASP;TOPMED=0.99998407237512742,0.00001592762487257
1   10851   rs1340958514    G   A   .   .   RS=1340958514;RSPOS=10851;dbSNPBuildID=151;SSR=0;SAO=0;VP=0x050000020005000002000100;GENEINFO=DDX11L1:100287102;WGT=1;VC=SNV;R5;ASP;TOPMED=0.99996018093781855,0.00003981906218144
1   10854   rs1222740568    G   A,C .   .   RS=1222740568;RSPOS=10854;dbSNPBuildID=151;SSR=0;SAO=0;VP=0x050000020005000002000100;GENEINFO=DDX11L1:100287102;WGT=1;VC=SNV;R5;ASP;TOPMED=0.99967348369011213,0.00019113149847094,0.00013538481141692
1   10857   rs1487378566    C   G   .   .   RS=1487378566;RSPOS=10857;dbSNPBuildID=151;SSR=0;SAO=0;VP=0x050000020005000002000100;GENEINFO=DDX11L1:100287102;WGT=1;VC=SNV;R5;ASP;TOPMED=0.99933900356778797,0.00066099643221202

However, it sounds like you might be underestimating the complexity here. There isn't just one minor allele, many sites have multiple minor alleles. The example above shows one such multiallelic site (you can split these into separate lines using bcftools). But you also have cases like:
1   10786   rs1464509424    CGGCGCA C   .   .   RS=1464509424;RSPOS=10787;dbSNPBuildID=151;SSR=0;SAO=0;VP=0x050000020005000002000200;GENEINFO=DDX11L1:100287102;WGT=1;VC=DIV;R5;ASP;TOPMED=0.99991239806320081,0.00008760193679918
1   10786   rs1258377526    CGGCGCAGGCGCAGACACATGCTAGCGCGTCGGGGTGGAGGCGT    C   .   .   RS=1258377526;RSPOS=10787;dbSNPBuildID=151;SSR=0;SAO=0;VP=0x050000020005000002000200;GENEINFO=DDX11L1:100287102;WGT=1;VC=DIV;R5;ASP

Don't expect just one major and one minor allele per location.

This is probably it, thanks. I looked on the dbSNP website but couldn't find the exact right file. I have a good method in place for handing multiallelic sites

I agree with @terdon that perhaps dbSNP is what you may want to look at. The dbSNP155 RefSNP VCF file in GRCh37 version may help. It can be downloaded from https://ftp.ncbi.nlm.nih.gov/snp/latest_release/vcf/GCF_000001405.25.gz.
You may also fine the 1000G genome website helpful. https://www.internationalgenome.org/data-portal/data-collection contains much of its release both in .fa format and in VCF format. I understand that a large number of samples included increase the file size tremendously, but it is necessary as you're not just looking at what the major and minor allele is, but also their frequency, which may vary greatly between populations. It's important to check which population we're seeing when using their genomic information.
Hope my answer helps. Regards.


text:
Does phasing improve imputation quality?
If so, by about how much? I want to know the downsides of skipping the phasing step on human microarray and WGS data.

What method do you want to use - this likely makes a difference.

@user438383 Beagle/Eagle

@user438383 but I'd also be interested in other methods if someone happens to know those too

Yes. Phasing data leads to a better imputation accuracy. More specifically:

Improves allele matching: Correct phased data ensures that alleles are matched correctly between the reference panel and the target dataset.
Using phased data, the computational time of imputation is significantly reduced.
Having haplotypes (phased data) helps imputing variants that were not genotyped (missing variants/genotypes).

This paper based on the UK10K data has some actual numbers describing imputation accuracy of phased vs unphased data.


text:
SNP How to detect regions of erratic pairwise alignment?
I have clusters of DNA sequences that contain some mutations. All the sequences in a cluster should contain the same mutations.
As the mutations aren't known, those clusters of sequences were pairwise-aligned against the non-mutated reference with Minimap2. Then an interesting fact cropped-up: a few mutations close to each other will make the alignment go awry in that zone.

Zone with a few SNPs:
ref  GTATGCGGATACCG

sq1  GTATCGTGACCG
sq2  GTATCTCGTGACCG
sq3  GTGTCTATTCGTGACCG
sq4  GTATCTCGTGACCG

Independent pairwise alignments against the reference (I removed the insertions from the output):
ref  GTATGCGGATACCG

sq1  GTAT-CGG--ACCG
sq2  GTATCTCGTGACCG
sq3  GTATTCGG--ACCG
sq4  GTATCTCGTGACCG

So I made a MSA of this small cluster (adding the reference in it) with Muscle5; the difference is noticeable, as we can clearly "see" the mutations now:

MSA alignment that includes the reference:
ref  GTATG---CGGATACCG

sq1  GTAT-----CGTGACCG
sq2  GTATC---TCGTGACCG
sq3  GTGTCTATTCGTGACCG
sq4  GTATC---TCGTGACCG

Same MSA, but trimmed:
ref  GTATGCGGATACCG

sq1  GTAT--CGTGACCG
sq2  GTATCTCGTGACCG
sq3  GTGTCTCGTGACCG
sq4  GTATCTCGTGACCG

I would say calculating the MSA seems to be the way to go but the clusters sizes range from 2 to 1000 sequences, the sequences are 3000nt long, and there's 2,000,000 sequences in total; we cannot afford to do the MSA of each cluster on the full size of the sequences.

On the other hand, we could probably manage to do a MSA of a few small zones in the sequences.
So, given a cluster and the pairwise alignments of its sequences against the reference, how would you determine the zones that need to be re-aligned?
note: I'm using Muscle5 because there's a Python binding of it and the results seem good. I also tried FAMSA, which is very fast, but it didn't yield good alignments.

Are sequences 1-4 in your 3 example supposed to be the same sequence? It is not really clear to me what you are trying to do / which problem you are trying to solve. You seem to be dealing mainly with indels rather than SNPs? Its also not clear why the alignment is awry, it is probably the best alignment with the given penalty matrix..

@Niklas Thank you for the comment. Yes, they are the same sequences (I updated the examples to make it clearer). By "awry" I mean that I cannot  make a consensus using the pairwise alignments, as they are subjectively "wrong" (mostly because of the indels)

OK. Generally, all alignments are equally "correct" or "wrong", since they minimize penalty score given the penalty matrix, which you can edit as you prefer. You would expect a msa to yield different results from pairwise alignments since it is a different problem.. regarding your actual question: I assume you could identify "bad" regions from a long stretch by the presence of many gaps in that region. But, keep in mind that when you re-align short sequences you may get different results than for the long sequences.

Maybe something like pileup is what you are looking for? https://samtools.github.io/bcftools/howtos/variant-calling.html

@Niklas Thank you, bfctool will definitely be useful


text:
How can I convert my tree to the correct nexus format?
I would like to use the program BayesTraitV4. I have my phylogenetic tree in NEXUS format but the program does not accept it. I checked the example data and while they are also in NEXUS format, they seemed different. I will provide a part of my data and the example data respectively.
My tree:
#NEXUS
begin taxa;
    dimensions ntax=85;
    taxlabels
    LevantinaFC4779
    Assyriella_naegelei
    Ceu_LAG_2
    CeuMist2
    CeuMist1
    Ceu_PAN_2
    Ceu_PAN_1
    Ceu_PAN_3
    Ceu_PAN_4
    CcMath1
    CcNav1
    CcNav2
    Cc_Mal_1
    Cc_KalR2_2
    CcRod1
    CcRod2
    CcRod3
    CcRod4
    CcHr1
    CeST6_1
    CeST6_4
    CeST6_5
    CeAA7_1
    CeST6_3
    CeST6_2
    Cg_MelR7_1
    Cg_MelR3_1
    Cg_MelR6_2
    Cg_MelR7_2
    Cg_MelR6_1
    Cg_MelR7_3
    Cg_MelR5_1
    Cg_MelR4_2
    Cg_MelR1_1
    Cg_MelR2_1
    Cg_MelR4_1
    Cg_MelR7_4
    Cg_MelR5_2
    ChTM4_2
    ChTM5_1
    ChTM4_3
    ChTM4_1
    Ch_MEN_5
    Ch_MEN_13
    Ch_MEN_8
    Ch_LEV_3
    Ch_LEV_5
    Ch_LEV_4
    Ch_MEN_6
    Ch_MEN_10
    Ch_LEV_2
    Ch_MEN_4
    Ch_MEN_14
    Ch_LEV_1
    Ch_MEN_11
    Ch_MEN_3
    Ch_MEN_9
    Ci_MSP_3
    Ci_MSP_2
    Ci_MSP_1
    Ci_DIA_1
    Ci_DIA_2
    CiPLAN1
    CiPLAN3
    CiPLAN2
    Ci_MEN_1
    Ci_MEN_2
    Ci_MEN_12
    CpGravR6_1
    CpGravR3
    CpGravR4
    CpGio_2
    CpGravR5_2
    CpGravR2
    CpGravR7
    CpGio_1
    CpGio_3
    CpGravR1
    CpGravR6_2
    CpGravR5_1
    CpThe_1
    CpSka_1
    Cn_IG_2
    Cn_IG_3
    Cn_IG_1
;
end;

begin trees;
    tree TREE1 = [&R] ((LevantinaFC4779[&rate_range={7.888601174550756E-4,0.2727277039812744},height_95%_HPD={0.0,1.5711876244495215E-12},length_range={0.6843282088933,6.197383612374},height_median=6.004086117172847E-13,length_95%_HPD={1.080900382604,3.875114824373},height=6.614724410322848E-13,rate=0.049900260878735954,height_range={0.0,5.799805080641818E-12},rate_median=0.045084434131968686,length=2.411967820837161,length_median=2.3472043999879997,rate_95%_HPD={0.0108358382756707,0.10141567250363237}]:2.4119678208371664,Assyriella_naegelei[&rate_range={2.4004562386929353E-4,0.19480057878691465},height_95%_HPD={0.0,1.5711876244495215E-12},length_range={0.6843282088933,6.197383612374},height_median=6.004086117172847E-13,length_95%_HPD={1.080900382604,3.875114824373},height=6.614724410322848E-13,rate=0.03027613080757391,height_range={0.0,5.799805080641818E-

Example tree:
#NEXUS
Begin trees;
    translate
        1   Opossum,
        2   Diprotodontian,
        3   Sloth,
        4   Armadillo,
        5   Anteater,
        6   Hedgehog,
        7   Mole,
        8   Shrew,
        9   Tenrecid,
        10  Golden_Mole,
        11  Sirenian,
        12  Hyrax,
        13  Elephant,
        15  Lo_Ear_Ele_shrew,
        17  Sciurid,
        18  Mouse,
        19  Rat,
        20  Hystricid,
        21  Caviomorph,
        22  Rabbit,
        23  Pika,
        24  Flying_Lemur,
        25  Tree_Shrew,
        26  Strepsirrhine,
        28  Phyllostomid,
        29  Free_tailed_bat,
        30  False_vampire_bat,
        31  Flying_Fox,
        32  Rousette_Fruitbat,
        33  Whale,
        35  Hippo,
        36  Llama,
        37  Ruminant,
        38  Pig,
        39  Horse,
        40  Rhino,
        41  Tapir,
        42  Cat,
        43  Caniform,
        44  Pangolin;
        tree tree.164000.208997.409491 = ((1:0.07628,2:0.065151):0.1546885,(((15:0.133631,(9:0.157087,10:0.077841):0.007403):0.005959,(11:0.0406,(12:0.092769,13:0.052284):0.001099):0.01453):0.029813,((4:0.05807,(3:0.054688,5:0.062757):0.008238):0.029441,(((26:0.083487,(24:0.064103,25:0.110063):0.003501):0.000775,((22:0.058663,23:0.116805):0.045375,(17:0.0898,((18:0.039267,19:0.043896):0.125808,(20:0.083297,21:0.089139):0.067194):0.006543):0.011062):0.005229):0.00803,((7:0.090648,(6:0.150383,8:0.133021):0.006503):0.01198,((36:0.062932,(38:0.066525,(37:0.076473,(35:0.039222,33:0.031429):0.00357):0.008228):0.005968):0.026714,((44:0.081521,(42:0.046822,43:0.059428):0.024607):0.004257,((39:0.046001,(40:0.029049,41:0.029533):0.003965):0.022397,((28:0.093611,29:0.050095):0.013432,(30:0.061026,(31:0.018601,32:0.018701):0.036144):0.004989):0.015808):0.001216):0.000872):0.002868):0.009567):0.010713):0.00234):0.1546885);
        tree tree.165000.208997.199794 = ((1:0.07628,2:0.065151):0.1546885,(((15:0.133631,(9:0.157087,10:0.077841):0.007403):0.005959,(11:0.0406,(12:0.092769,13:0.052284):0.000485):0.01453):0.029813,((4:0.05807,(3:0.054688,5:0.062757):0.008238):0.029441,(((26:0.083487,(24:0.064103,25:0.110063):0.002826):0.000775,((22:0.058663,23:0.116805):0.045375,(17:0.0898,((18:0.039267,19:0.043896):0.125808,(20:0.083297,21:0.089139):0.067194):0.006543):0.011062):0.003736):0.00803,((7:0.090648,(6:0.150383,8:0.133021):0.010132):0.01198,((44:0.084486,(42:0.046822,43:0.059428):0.024607):0.004257,((36:0.062932,(38:0.066525,(37:0.076473,(35:0.039222,33:0.031429):0.004674):0.008228):0.005968):0.026714,((39:0.046001,(40:0.029049,41:0.029533):0.003965):0.022397,((28:0.093611,29:0.050095):0.008576,(30:0.061026,(31:0.018601,32:0.018701):0.036144):0.004989):0.015808):0.001216):0.000872):0.002868):0.009567):0.010713):0.00234):0.1546885);
;

end;

Why are these formats different and how can I make my file to look like the example one?
Keep in mind the files are not presented in their entirety because I'm limited with the amount of characters.
Edit: I will be providing my tree in this link https://drive.google.com/file/d/1lYmzW9htwaclC_yYs3fQom-Z3zbjk8wH/view?usp=sharing

Okay ... this will probably work
#NEXUS
begin trees;
    tree TREE1 = [&R] ((LevantinaFC4779:2.411968,Assyriella_naegelei:2.411968):5.752052,((Cn_IG_1:0.126901,(Cn_IG_2:0.075971,Cn_IG_3:0.075971):0.05093):7.4805,(((CpThe_1:0.122189,CpSka_1:0.122189):2.663216,((CpGio_2:0.137489,(CpGravR6_1:0.087156,(CpGravR3:0.038234,CpGravR4:0.038234):0.048922):0.050334):0.257308,(CpGravR5_1:0.28159,((CpGravR7:0.099712,(CpGravR5_2:0.020255,CpGravR2:0.020255):0.079458):0.15617,(CpGravR6_2:0.109222,(CpGio_1:0.068711,(CpGio_3:0.049345,CpGravR1:0.049345):0.019367):0.040511):0.14666):0.025708):0.113207):2.390607):1.772627,(((Ci_MEN_12:0.072368,(Ci_MEN_1:0.022725,Ci_MEN_2:0.022725):0.049643):1.78894,((CiPLAN2:0.253487,(CiPLAN1:0.067509,CiPLAN3:0.067509):0.185978):0.386279,((Ci_MSP_3:0.061614,Ci_MSP_2:0.061614):0.231268,(Ci_MSP_1:0.113777,(Ci_DIA_1:0.047073,Ci_DIA_2:0.047073):0.066704):0.179105):0.346885):1.221542):1.580745,(((Ceu_LAG_2:1.580984,(CeuMist2:0.198982,CeuMist1:0.198982):1.382002):0.341226,((Ceu_PAN_2:0.046898,Ceu_PAN_1:0.046898):0.20368,(Ceu_PAN_3:0.018814,Ceu_PAN_4:0.018814):0.231764):1.671632):1.226475,((CcHr1:1.247104,((CcRod1:0.403751,(CcRod4:0.276104,(CcRod2:0.158977,CcRod3:0.158977):0.117127):0.127647):0.308704,((Cc_Mal_1:0.138232,Cc_KalR2_2:0.138232):0.395435,(CcMath1:0.338389,(CcNav1:0.091762,CcNav2:0.091762):0.246627):0.195278):0.178788):0.534649):1.609976,(((CeST6_5:0.056918,(CeST6_1:0.030107,CeST6_4:0.030107):0.026812):0.733078,((CeST6_2:0.143014,(CeAA7_1:0.097814,CeST6_3:0.097814):0.0452):0.349431,((Cg_MelR7_1:0.136715,(Cg_MelR3_1:0.088065,(Cg_MelR6_2:0.047987,(Cg_MelR7_2:0.019268,Cg_MelR6_1:0.019268):0.028719):0.040078):0.04865):0.074941,(Cg_MelR5_2:0.143401,(Cg_MelR7_3:0.119131,(Cg_MelR5_1:0.090065,((Cg_MelR4_1:0.020938,Cg_MelR7_4:0.020938):0.051613,(Cg_MelR4_2:0.052913,(Cg_MelR1_1:0.021387,Cg_MelR2_1:0.021387):0.031526):0.019639):0.017514):0.029066):0.02427):0.068255):0.280789):0.297551):0.793423,(((ChTM4_2:0.049871,ChTM5_1:0.049871):0.111934,(ChTM4_3:0.023134,ChTM4_1:0.023134):0.138672):0.228306,(Ch_MEN_9:0.174548,((Ch_MEN_11:0.021817,Ch_MEN_3:0.021817):0.10954,((Ch_MEN_8:0.057141,(Ch_MEN_5:0.022973,Ch_MEN_13:0.022973):0.034168):0.053065,((Ch_LEV_3:0.031056,Ch_LEV_5:0.031056):0.061225,((Ch_LEV_4:0.029359,(Ch_MEN_6:0.017456,Ch_MEN_10:0.017456):0.011904):0.034407,((Ch_LEV_2:0.017911,Ch_MEN_4:0.017911):0.040741,(Ch_MEN_14:0.03419,Ch_LEV_1:0.03419):0.024462):0.005114):0.028515):0.017925):0.021151):0.043191):0.215563):1.193307):1.27366):0.291605):0.293368):1.115978):3.04937):0.556619);
end;

Alternatively, I'm not sure about the root, just remove [&R] in the above.
Finally there could be an issue that you haven't obtained the sample directly from the raw Beast output, instead you've summarised it. It might complain about that. Thus it might need all those trees from the raw Beast output. Thats not a format issue it's a data input issue. BayesTraits might be perform the calculation against the total MCMCMC sample minus the burn-in.

from comments
Read it into Rs ape::read.nexus() and then write via phytools::writeNexus(). This should produce your translate table. Biopython recognises this format but doesn't write to it.
The closest Biopython can achieve is a taxablock here
#NEXUS
Begin Taxa;
 Dimensions NTax=85;
 TaxLabels LevantinaFC4779 Assyriella_naegelei Cn_IG_1 Cn_IG_2 Cn_IG_3 CpThe_1 CpSka_1 CpGio_2 CpGravR6_1 CpGravR3 CpGravR4 CpGravR5_1 CpGravR7 CpGravR5_2 CpGravR2 CpGravR6_2 CpGio_1 CpGio_3 CpGravR1 Ci_MEN_12 Ci_MEN_1 Ci_MEN_2 CiPLAN2 CiPLAN1 CiPLAN3 Ci_MSP_3 Ci_MSP_2 Ci_MSP_1 Ci_DIA_1 Ci_DIA_2 Ceu_LAG_2 CeuMist2 CeuMist1 Ceu_PAN_2 Ceu_PAN_1 Ceu_PAN_3 Ceu_PAN_4 CcHr1 CcRod1 CcRod4 CcRod2 CcRod3 Cc_Mal_1 Cc_KalR2_2 CcMath1 CcNav1 CcNav2 CeST6_5 CeST6_1 CeST6_4 CeST6_2 CeAA7_1 CeST6_3 Cg_MelR7_1 Cg_MelR3_1 Cg_MelR6_2 Cg_MelR7_2 Cg_MelR6_1 Cg_MelR5_2 Cg_MelR7_3 Cg_MelR5_1 Cg_MelR4_1 Cg_MelR7_4 Cg_MelR4_2 Cg_MelR1_1 Cg_MelR2_1 ChTM4_2 ChTM5_1 ChTM4_3 ChTM4_1 Ch_MEN_9 Ch_MEN_11 Ch_MEN_3 Ch_MEN_8 Ch_MEN_5 Ch_MEN_13 Ch_LEV_3 Ch_LEV_5 Ch_LEV_4 Ch_MEN_6 Ch_MEN_10 Ch_LEV_2 Ch_MEN_4 Ch_MEN_14 Ch_LEV_1;
End;
Begin Trees;
 Tree tree1=((LevantinaFC4779:2.41197,Assyriella_naegelei:2.41197):5.75205,((Cn_IG_1:0.12690,(Cn_IG_2:0.07597,Cn_IG_3:0.07597):0.05093):7.48050,(((CpThe_1:0.12219,CpSka_1:0.12219):2.66322,((CpGio_2:0.13749,(CpGravR6_1:0.08716,(CpGravR3:0.03823,CpGravR4:0.03823):0.04892):0.05033):0.25731,(CpGravR5_1:0.28159,((CpGravR7:0.09971,(CpGravR5_2:0.02025,CpGravR2:0.02025):0.07946):0.15617,(CpGravR6_2:0.10922,(CpGio_1:0.06871,(CpGio_3:0.04934,CpGravR1:0.04934):0.01937):0.04051):0.14666):0.02571):0.11321):2.39061):1.77263,(((Ci_MEN_12:0.07237,(Ci_MEN_1:0.02272,Ci_MEN_2:0.02272):0.04964):1.78894,((CiPLAN2:0.25349,(CiPLAN1:0.06751,CiPLAN3:0.06751):0.18598):0.38628,((Ci_MSP_3:0.06161,Ci_MSP_2:0.06161):0.23127,(Ci_MSP_1:0.11378,(Ci_DIA_1:0.04707,Ci_DIA_2:0.04707):0.06670):0.17910):0.34688):1.22154):1.58075,(((Ceu_LAG_2:1.58098,(CeuMist2:0.19898,CeuMist1:0.19898):1.38200):0.34123,((Ceu_PAN_2:0.04690,Ceu_PAN_1:0.04690):0.20368,(Ceu_PAN_3:0.01881,Ceu_PAN_4:0.01881):0.23176):1.67163):1.22647,((CcHr1:1.24710,((CcRod1:0.40375,(CcRod4:0.27610,(CcRod2:0.15898,CcRod3:0.15898):0.11713):0.12765):0.30870,((Cc_Mal_1:0.13823,Cc_KalR2_2:0.13823):0.39543,(CcMath1:0.33839,(CcNav1:0.09176,CcNav2:0.09176):0.24663):0.19528):0.17879):0.53465):1.60998,(((CeST6_5:0.05692,(CeST6_1:0.03011,CeST6_4:0.03011):0.02681):0.73308,((CeST6_2:0.14301,(CeAA7_1:0.09781,CeST6_3:0.09781):0.04520):0.34943,((Cg_MelR7_1:0.13672,(Cg_MelR3_1:0.08807,(Cg_MelR6_2:0.04799,(Cg_MelR7_2:0.01927,Cg_MelR6_1:0.01927):0.02872):0.04008):0.04865):0.07494,(Cg_MelR5_2:0.14340,(Cg_MelR7_3:0.11913,(Cg_MelR5_1:0.09007,((Cg_MelR4_1:0.02094,Cg_MelR7_4:0.02094):0.05161,(Cg_MelR4_2:0.05291,(Cg_MelR1_1:0.02139,Cg_MelR2_1:0.02139):0.03153):0.01964):0.01751):0.02907):0.02427):0.06825):0.28079):0.29755):0.79342,(((ChTM4_2:0.04987,ChTM5_1:0.04987):0.11193,(ChTM4_3:0.02313,ChTM4_1:0.02313):0.13867):0.22831,(Ch_MEN_9:0.17455,((Ch_MEN_11:0.02182,Ch_MEN_3:0.02182):0.10954,((Ch_MEN_8:0.05714,(Ch_MEN_5:0.02297,Ch_MEN_13:0.02297):0.03417):0.05307,((Ch_LEV_3:0.03106,Ch_LEV_5:0.03106):0.06123,((Ch_LEV_4:0.02936,(Ch_MEN_6:0.01746,Ch_MEN_10:0.01746):0.01190):0.03441,((Ch_LEV_2:0.01791,Ch_MEN_4:0.01791):0.04074,(Ch_MEN_14:0.03419,Ch_LEV_1:0.03419):0.02446):0.00511):0.02851):0.01792):0.02115):0.04319):0.21556):1.19331):
1.27366):0.29161):0.29337):1.11598):3.04937):0.55662):0.00000;
End;

The solution could be coded from the taxa block if the above doesn't work, but easier to use R ape and phytools

I tried using this tree you provided as input and got this message:
Err: Does have a translate line.


text:
INDELS in PLINK files converted to VCF
I want to compare/validate variants called from sequencing data with array (plink format) variant data. I converted the plink files (.bim, .bed, and .fam files) with plink1 to vcf files.
plink --bfile prefix_plink --recode vcf-iid --out prefix_out

However, the plink vcf files have "I" and "D" values for INDEL variants that cannot be compared and break the tools downstream (I am now using hap.py with rtgtools to validate). Is there a way to convert plink files to vcf, with some transformation of these I and D variants to a usable value? A "D" I could translate by hand to "" which is a vcf accepted format. But insertions require I guess a more complex procedure, or are even impossible (because they need to be compared to a reference genome)? Or am I missing some flag somewhere?
I would love to get some pointers and tips for this.
Kind regards,

I believe bcftools is able to accomplish this.


text:
Is 16S sequencing cheaper than whole-genome sequencing for taxonomic assignment?
This is adapted from a question posted on reddit
In my work, I occasionally get asked about metagenomic sequencing, which (after further questioning) I discover means that the person wants to work out the taxonomic composition of a sample.
In the past, people have gone straight to 16S sequencing for this, and various manufacturers have created targeted 16S kits for this purpose. But does it still make sense to do that? Is 16S still cheaper and better than other options?

When comparing anything targeted to whole-genome sequencing for the purposes of assembly, the targeted approach will always be cheaper, because the total length that needs to be assembled is shorter. But that's like comparing a video to a photograph.
If you're just looking to determine the phylogenetic composition of a bacterial sample, then there's no reason why targeted 16S sequencing should be cheaper than whole-genome shotgun sequencing, and the whole genome approach may end up better due to the highly conserved nature of the 16S sequence, and the abundance of horizontal gene transfer.
Previously, 16S might have been better because the taxonomic databases were more complete, but now that there are so many microbial genomes available, the whole genome databases are comprehensive enough for most purposes.
@terdon linked to this 2008 paper, in which they compare random sequence reads to 16S for estimating the biodiversity of a metagenomic library, finding that 16S wasn't worth it and they could get perfectly good results with random reads.
In a situation where samples are abundant and no amplification is required, PCR can be excluded from whole-genome sequencing, making it the cheaper option at the sample preparation side of things.
Just to put numbers to things to demonstrate how cheap whole-genome sequencing can be, here are the big consumable costs for nanopore (prices in USD):

RBK114.24 rapid barcoding kit* - \$700 (24 samples, 12 runs)
Flongle flow cell - \$810 (12 flow cells)

So the marginal cost of one run is \$68 for the flow cell, and \$60 for the sequencing kit, processing 24 samples at a cost of about $11 per sample (rounding up).
[FWIW, ONT also provide a rapid PCR barcoding kit for the same cost (i.e. \$700, but polymerase is not included), which can be used together with PCR for low-yield samples.]
Flongle flow cells will usually generate >200 Mb per run (up to about 2 Gb in rare cases), so that's about 8 Mb of sequence per sample, which should be sufficient to do proportional counting down to about 1% of the sample (assuming >1000 reads per sample).
* There is a 96-barcode kit available for \$990, but I don't think that would have enough coverage to be useful when running samples on a Flongle flow cell.


text:
What does the different sequences represent?
I am using this package nsdpy to download genome sequences from NCBI nucleotide database.
Specifically I am interested in the whole mitochondrial genome of different species, here I will use a subset of species just as an example.
My Query is:
"{mitochondrion[Title]" -i -t -L subset_species_list.txt
Where :

-i Adds the taxonomic information to the information lines of the sequences written in the output files (not relevant for this example)
-t option produces a tsv file with the following fields: Taxon name, SequenceID, TaxID, Taxonomic lineage, Sequence length, Sequence.
subset_species_list.txt is a list of species, below you can see its content

Xenopus laevis
Gasterosteus aculeatus
Chironomus riparius
Rasbora heteromorpha
Planktolyngbya limnetica

The output file is quite long and includes the actual sequences. Here I will post the same file where I dropped the column for the sequence and the lineage to keep it compact

My aim is to obtain phylogenetic distances between species.
My questions are:

What do the sequences for the same species represent?
Which sequences are best suited for an alignment and how do I know that?
Should I also look into the orientation of the sequences or are they all on the forward strand?
Is it better if I filter for Coding sequences (CDS) instead of whole genomes

Also if you noticed any rookie mistake in my process feel free to point that out. It is the first time for me working with DNA sequences and I really want to learn to correct way to do that.

Welcome to the site. This question is too broad could you limit the number of questions please?

Welcome! I have some quick basic answers for the questions you're asking so hopefully this helps:

Sequences for the same species represent variation/diversity/evolutionary changes. Think of humans; how do different people's DNA sequences represent the sequence of a human? This is especially important for phylogenetic studies.
For the most part, all sequences are best suited for alignment. The main thing to think about is the amount of data you have (if it's a lot of GB the alignment will take a lot of computing power), the type of sequnece data you have (simple rule of thumb is to not mix short read with long read sequences), and the tool you're using for alignment.
Always a good idea to think about orientation but most tools will account for orientation when aligning sequences.
This depends on what you're looking for but more often than not using the whole genomes is best as it lets you look at the whole picture instead of just a small portion but that's definitely a lot of computing power. If you're just looking at the mitochondria just focus on that region

Hope this helps!