text:
Alignment with arbitrary number of mismatches or gaps
I have 23bp long reads and want to find all possible alignments of them to the human genome (hg19, hg38) for an arbitrary number of mismatches (<7), possibly also small indels. I've read in literature that people use bowtie2 for this, so I've tried the following:
echo "GGAGAACAAGTGGCAGATAGAGG" | bowtie2 --all -f -U- -x hg19 -r --no-hd --quiet --score-min C,-36 -p 3 -L 4

I'm not very experienced with bowtie2, and I hoped that this would give me all alignments of the piped-in sequence with up to 6 mismatches (6*-6 = -36) to the hg19 genome which I previously indexed from fasta files. It however only gives me the correct full alignment:
0   0   chr2    173327389   255 23M *   0   0   GGAGAACAAGTGGCAGATAGAGG IIIIIIIIIIIIIIIIIIIIIII AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:23 YT:Z:UU

I've read that bowtie2 by design only allows a maximum of one mismatch per seed. Even if I implement the smallest possible seed length and interval (-L 4 -i C,1 -N 1), I only get this one result which I think is highly unlikely, given that I'm aligning to the whole genome.
The last three bp will have to match exactly - I was going to enforce this manually later on, maybe there's a better option for doing this. Taking the indels into account is also a secondary aim.
I'm fine with ~10k alignments per query for <=6 mismatches plus the last 3bp matching exactly. Haven't done the combinatorics but that sounds realistic to me - question is how to make bowtie do it, or a different tool.
What would be the right way of achieving this using bowtie, or a different tool?

This is not something mappers will be good at, because you're going to hit most of the genome. Bowtie won't do it, and I don't think LAST will either. The majority of fast mapping programs rely on seed-based mapping where the seed is essentially free of errors.
Doing only mismatches makes it much easier to determine, but it's still too relaxed. Consider that the example sequence you have provided is 10 mismatches away from GGGGGGGGGGGGGGGGAGAGAGG, and also 10 mismatches away from AAAAAAAAAAAAAAAGAGAGAGG (with the last 3 bp matching exactly)
Depending on how many reads you have it may be better to work in reverse, iterating over all the kmers in the genome and working out whether or not they are fewer than 7 mismatches (plus indels) away from any reads:
for gp in genomePositions
    genomeKmer = genome[gp:(gp+23-1)]
    for read in reads
        errors = 0
        for rp in 1..23
            if (genomeKmer[rp] != read[rp])
                errors += 1
        if (errors < 7)
            print gp

I've ended up using a package called BatMis for this, which can detect alignments with an arbitrary number of mismatches (not sure about indels). It is quite performant and yields about the order of magnitude of alignments that I was expecting (several thousands for alignment of a 20bp read to the whole human genome with up to 6 mismatches).

Did you see

-k  Report up to  valid alignments per read or pair (default: 1).

Did see that, but I thought --all (or -a in short) means it outputs all found alignments, ie supersedes -k? Using -k 5 instead of --all still only yields one result.


text:
maf file format: Variant_Classification value for non frameshift substitution
I have a vcf containing SNV and MNV (short indells or substitutions), I have annotated that file with annovar and in the resulting txt file I have this row:
chr7    140453136   140453137     AC   TT  exonic   BRAF    .     nonframeshift substitution              BRAF:NM_004333:exon15:c.1798_1799AA:p.V600K

If I convert the txt file in maf using annovarToMaf function of maftools bioconductor package I obtain NA as value for Variant_Classification field. In other words, the field is included in the original vcf file but it is converted into NA when saved as maf with maftools.
What is the correct value of the Variant_Classification field in a .maf file for non-frameshift substitutions like BRAF V600K?
I tried vcf2maf and that tool put Missense_Mutation as a value for Variant_Classification of BRAF V600K mutations.
I tried to download the .maf file of TCGA melanoma (that should contain BRAF V600K mutations) but the access to that file is restricted.

Ok, try to make a reproducible example as it might be a bug on the maftools package. However I would first post this reproducible example here to see if there is a solution without involving changes on the package

source: Stack Exchange

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

sort: votes


text:
Improving list speed: glycosylation example
I discussed a question with @gaspanic Python/Biopython - Replace amino acid residue on MSA with "z" from a list of unaligned positions
. The issue emerged was speeding up lists in Python. The solutions we discussed were

Raw code
Predefining a list size
deque

1. @gaspanic
my code
2. @gaspanic
my code
3. @gaspanic
my code

3.45-3.48
7.17-7.2
4.38-4.85
7.09-7.23
3.44-3.53
7.21-7.39

Run timings are seconds (3 replicated) above with @gaspanic's code being faster.
What improvements would speed this up?
Here I extended the data set
from itertools import repeat
def iterit(x):
    seq = 'M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K'
    return ''.join(list(repeat(seq, x)))

if __name__ == '__main__':
    seqE3 = iterit(2800) # length is 201600
    seqE4 = iterit(28000) # length is 2016000
    seqE5 = iterit(280000) # length is 20160000

@gaspanic code
for seqbit in [seqE3,seqE4,seqE5]:
    seq = seqbit
    pos = [2, 5, 10, 11, 14]
    new_seq = []
    count = 0
    for s in seq:
        if s == "-":
            new_seq.append(s)
        else:
            count += 1
            if count in pos:
                new_seq.append('Z')
            else:
                new_seq.append(s)
   # print(''.join(new_seq))

My code:
# same variables and loop structure
newseq = list()
count = 0
for n,s in enumerate(prot):
    if str(s) == '-':
        count += 1
        s = '-'
    else:
        if (n + 1 - count) in pos:
            s = 'Z'
    newseq.append(s)
# print(''.join(newseq))

Predefining lists code (additional counter needed for @gaspanic's code:
# same variables and loop structure
new_seq = [] * len(seq)
count = 0
for n,s in enumerate(prot):
    if str(s) == '-':
        count += 1
        s = '-'
    else:
        if (n + 1 - count) in pos:
            s = 'Z'
    newseq[n] = s

Deque code simply any of the above substituting
newseq = list()

with
from collections import deque
newseq = deque()

Suggestions for efficiency improvements welcome.
A solution in Cython code would be cool. PyPI is not sought here.

I'd start by saying that Python is inherently slow and this feels like a case of perhaps unnecessary optimisation. With that in mind here are my suggested (untested) optimisations of @gaspanic code:
pos = {2, 5, 10, 11, 14} # make pos a set for O(1) member testing, and don't define inside the loop
for seq in (seqE3,seqE4,seqE5): # Define seq here rather than inside the loop, and use an immutable tuple instead of a mutable list
    new_seq = []
    count = 0
    for s in seq:
        if s == "-":
            new_seq.append(s)
        else:
            count += 1
            new_seq.append('Z' if count in pos else s) # use ternary operator

Thanks, the tuple did add speed, particularly tuple the iterator output. Using a set (the {} brackets) was a good solution and reduced runtime by ~0.5 seconds, which I'd not considered. It was better than tuple in this example. Repositioning pos adds a tiny amount, but agree its better practice and depends on the size of pos (in this case its small). Happy to discuss Python's slow speed in a summary.

Huge thanks to @Chris_Rands, this solution increases runtime by ~ 0.5 seconds uniformly across all coded solutions, so this code has a runtime of just over 3 seconds runtime on the benchmark machine. I incorporated my solution with @Chris_Rands solution and the final runtime I got was:
@gaspanic

1.7 - 1.8 seconds

This being twice the fastest speed in the question and four times the speed of my raw code.
My code

4.4 - 4.5 seconds

However a unique combination is needed:

tuple (Chris Rands)
set - important (Chris Rands)
an iterator
a list comprehension

In particular 3. and 4. are needed in combination. This is explained below.
def iterit(x):
    seq = 'M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K'
    return tuple(''.join(list(repeat(seq, x))))

def mIterit(seqE3,seqE4,seqE5):
    pos = {2, 5, 10, 11, 14}
    for seqbit in (seqE3,seqE4,seqE5):
        seq = seqbit
        count = 0
        # loc = 0
        for s in seq:
            if s == "-":
                yield '-'
            else:
                count += 1
                if count in pos:
                    yield 'Z'
                else:
                    yield s

f __name__ == '__main__':

    seqE3 = iterit(2800) # length is 201600
    seqE4 = iterit(28000) # length is 2016000
    seqE5 = iterit(280000) # length is 20160000

    listComp = [i for i in mIterit(seqE3,seqE4,seqE5)]

The list comprehension (variable listComp) is the very last line and is crucial to the solution otherwise there is little increase in speed to @Chris_Rand's code. This alternative is to replace the last line with:
newseq = list()
for i in mIterit(seqE3,seqE4,seqE5):
    newseq.append(i)

A simple:
for i in mIterit(seqE3,seqE4,seqE5):
    pass

Is the same speed as the list comprehension, but its of no value because the output isn't saved to anything.
I was impressed by the interaction between the iterator and list comprehension. The other notable influence on speed is the set described by @Chris_Rand.
How does it work? tuples do increase speed by removing some functional of a list (immutability). Here the gain is small. Placing the pos variable outside the loop speaks for itself (0.1 seconds gain). set is interesting.
The iterator spews out a stream of results via yield which is where it attains its speed. The problem is capturing that stream outside the iterator is expensive. However list comprehensions are faster than standard loops and in this scenerio carry little/no detectable overhead. I was impressed. Up to this point I thought list comprehensions were a bit faster and a lot 'cooler' than loops. I hadn't understood there was a synergy with an iterator.
Is it all worth it? Python can easily achieve faster run times than Java using Cython but I will place the discussion in the chat forum next week.

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


text:
Refactoring pandas using an iterator via chunksize
This question was also asked on Stack Overflow
Bioinformatics rationale eggNOG files can be very big and sump all available RAM for regular to medium sized desktops.

I am looking for advice on using a pandas iterator.
I performed a parsing operation using Python pandas, the size of the input files is  eggNOG a program used in metagenomic gene function identification but has broader application. The files are large and parsing is resulting in 'RAM bottleneck' phenomenon and will not processing/parse the file.
The obvious solution is to shift to an iterator, which for pandas is the  chunksize option
import pandas as pd
import numpy as np

df = pd.read_csv('myinfile.csv', sep="\t", chunksize=100)

Whats changed with the original code is the chunksize=100 bit, forcing an iterator.
The next step is just to perform a simple operation, dropping a few columns and moving all '-' characters to np.nan then writing the whole file to disk.
df.drop(['score', 'evalue', 'Description', 'EC', 'PFAMs'],axis=1).replace('-', np.nan)
df.to_csv('my.csv',sep='\t',index=False)

How is the step performed under an iterator?

The solution posted by @el_oso and @Corralien here was as follows:

if you're already having memory issues reading the file, the below will actually read and write to the new file in chunks.

cols_to_keep = {'col1', 'col2', 'col3'}
df = pd.read_csv('myinfile.csv', chunksize=100, usecols=cols_to_keep)
with open('my.csv', 'w') as f:
    f.write('\t'.join(cols_to_keep) + '\n')
    for chunk in df:
        #do your processing here, appending that chunk to the file
        chunk.to_csv(f, header=None, sep='\t')

f is _io.TextIOWrapper i.e. a stream you can consistently write to in small memory chunks

This is very elegant because there are two components which are both cool and highlight why pandas is so widely used.

Instead of loading a large file and then dumping stuff thats not wanted: don't load junk at point source.
Simply open the write file outside the chunking loop will then stream the chunks to straight to disk and once the final chunking loop completes the file snaps shut.

The process is

Load a chunk of data, but don't load any 'junk'.

Process the chunk of data and write it to disk.

Load another chunk of data repeating the above two steps.

Continue until the entire file is processed.

The "do processing" is an OO module (there's alot going on) and the chunk is actually self.chunk. None of the original code needs to change its simply restructured inside the open loop to stream chunks to disk.
The beauty of this solution is chunksize is a variable, here 100 is just the example. Python's iterator called yield is one line at a time iteration, but this the automates a defined 'chunk' at a time. Thus the runtime speed can be tuned against the "RAM bottleneck".
Testing The chunked outfile was then reduced via groupby to remove the "duplicates" and frequencies summed (code not shown). The chunking OO module output file is identical to "non-chunked" OO module and written back into the module.
Its worth upvoting the SO answers.


text:
How to visualize electron density of a specific part of the protein residue?
Hi everyone (mostly who are interested in structural bioinformatics).
I face the problem that I cannot plot the electron density of a specific part of my protein. For the whole protein it is not a problem:
from Bio.PDB import *
import nglview as nv
import ipywidgets

pdbl = PDBList()
path = pdbl.retrieve_pdb_file('1FAT')
parser = MMCIFParser()
struct = parser.get_structure("1fat", path)
view = nv.show_biopython(struct)
view

But when I tried to pick up residues from 10 till 109 (res_list = Selection.unfold_entities(struct, "R")[10:109]) and plot I can't do that, because for nglview I need to have a special Object. As I'm a newbie, I'd like to address this question to specialists.
At the moment I'm trying to visualise model (structure), not electron density. But the next step is trying to find a package to visualize electron density also. I've just started becoming familiar with some open packages.
Thank You in Advance!

Hi @KitsuneBreeze its a good question and looks to be resolved. Please remember to "accept" and "upvote"

NGLview runs off the JavaScript libray NGL, which is documented as a series of very helpful examples: in the gallery. The commands are not one-to-one (slightly frustratingly), but it is helpful as you can figure it out easily with a few dir and help commands. In this case whatever example of density being loaded it will reveals that it's a different component —this is the same everywhere as density files do not contain atomic coordinates. Unfortunately, there's not a direct way to generate a new ComponentViewer from different sources to an NGLWidget and you cannot merge two of them: the various .show_* functions return an NGLWidget but aren't part of it.
First, I am not sure why you are going through a convoluted route for loading a PDB code:
pdb_code = '1fat'

import nglview as nv

view: nv.NGLWidget = nv.show_pdbid(pdb_code)

Now, on loading the component a representation is added to it.
These can hide the representation you then add and things will not make sense, until this detail is recalled! I fall in that pit nearly every time.
But we can change it with update command:
# I don't like wonky sheets
view.update_cartoon(smoothSheets=True, color='salmon', component=0)

Now let's add the density map the circuitous way (no .show_url).
import requests, io
response: requests.Response = requests.get('https://www.ebi.ac.uk/pdbe/coordinates/files/1fat.ccp4')
response.raise_for_status()
fh = io.BytesIO(response.content)
comp: nv.component.ComponentViewer = view.add_component(fh, ext='ccp4')

Now we can add a surface representation, but first let's clear the default.
For the arguments I had a gander at the example:
view.clear_representations(component=1)
view.add_representation('surface',
                        color='grey',
                        wireframe=True,
                        isolevel=1.5,
                        component=1)
view.center('100-110:A')
view


text:
Extracting Snakemake metadata (filenames) with the Snakemake API
I am creating a program that can store all relevant data from a Snakemake analysis run into a gzipped file (or a Pickle, not sure yet), for easy backup and moving around different machines.
To do this, it would be very helpful I I could somehow use my samples.tsv file and the Snakefile (the rules) to create a list of all files that are created in the entire data processing pipeline. I just want a list of all files created, per sample, during a workflow.
I was wondering if that is perhaps possible to do with the Snakemake API. I see some functions that seem useful but I don't quite understand the documentation without examples.
Can anyone help? Are there perhaps other ways to achieve the same thing? I could write something myself of course (parse all filenames from the rules), but I think this must have been created before...
Thank you.
Edit: Providing some more details of what I am looking for, for example I would like to do this:
import snakemake.api

workflow = snakemake.api.SnakemakeApi.workflow(
  'path/to/Snakefile',
  'path/to/config',
  'path/to/samples.tsv'
  )

workflow.files_created_in_workflow

>

[
  'expression/SRR29098682.genes.results',
  'expression/SRR29098675.genes.results',
  'expression/SRR29098683.genes.results',
  'qc/SRR29098682.rseqc.fragment_length.txt'
 ...
]

Etc, so that I can then use that list of files to gzip them all together (after excluding some large files like .bam and .fastq.gz), preferrably on a per sample basis.


text:
BAM to gene expression matrix (UMI counts per gene per cell),10X
I am trying to reproduce some results of a scRNASeq experiment. However I am new to the server-side aspect of such analyses and am very confused at the moment.
The data provided by the authors of the paper is in .BAM format and from there I wish to derive a gene expression matrix (UMI counts per gene per cell). The authors stated that they did this using the 10X genomics pipeline, however it seems that .BAM outputs are sort of the "final version" of the pipeline and, the pipeline only really deals with truly raw sequencing files (bcl/fastq). I am trying to learn how to do all this processing myself as I will need to eventually. Unfortunately, most of the documents I have read end up making me more confused. I know BAM files are technically already processed and aligned but, how would I get the counts from the file?
I have considered converting it back to fastq format and following the pipeline from the beginning, but it just seems like I would be moving backwards by doing such.
I have also found some R packages that can integrate BAM files in the environment like Rsamtools, but I would much rather do this on the server and then upload a csv into R for the downstream analysis as the files are +60GB each.
Essentially I am just asking for some advice, or a vignette/article I can read that would clear up how to proceed with the BAM outputs to derive the gene expression data. I have already read the cellranger pages on the 10X genomics website but they were of no help in answering this quesion.
Thank you for your time,

How to re-analyze 10X BAM files?
This is a great question and honestly, I don't think there was an easy way to do this at the time the question was asked. The reason is that if you want to re-do the authors' analysis to get the gene expression (gene-cell / feature-barcode) matrix and they used cellranger (the official 10X pipeline software) to do that, you will need FASTQ files as input - because "10x pipelines require sequencer FASTQs (with embedded barcodes) as input."
Therefore, to get FASTQ files from 10X BAM files, there's a relatively new tool provided by 10X Genomics: bamtofastq. Then you can simply use those FASTQ files as input to cellranger.
More info here: https://support.10xgenomics.com/docs/bamtofastq
How to correctly download 10X BAM files from NCBI?
The only caveat is that you need a BAM file generated directly by 10X's cellranger (or the respective 10X pipeline, if not dealing with gene expression) - that means that a BAM file obtained by downloading an SRA from NCBI and converting to BAM won't work; you need to get the original BAM file directly (often found among the originally submitted files, under "Original format").
To give you some additional help here (because figuring this out is also not trivial): an original 10X BAM file (if originally submitted in this format) can be downloaded from NCBI using prefetch from sra-tools, specifying the --type option - for example:
prefetch --type TenX --max-size 100000000 SRR5167880

Because the file was larger than the default limit of 20GB, I also had to increase it by specifying the --max-size option.
Alternatively, you can follow the official 10X Instructions to Download and Process BAM files of 1.3 Million Brain Cells.
I know this is an old thread but I'm hoping my answer will be helpful to anyone who has the same question (like me).

This is the simplest approach (if the process’s data is not available on GEO). Computing the UMI count matrix for 10x Genomics data is already implemented in Cell Ranger and if the same version is used it should reproduce the results of the study.

You should be able to parse out what you need using the tags in the .bam.  10xGenomics' website says what tags they add.
https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/bam
Going backwards would also involve parsing the tags, because you have to make two fastq files, and a simple bam -> fastq pipeline won't do that correctly.

Could you be more specific please? I'm very new to this procedure and bash scripting in general.

Yes, 10xGenomics single cell bams have the genes called in the bam tags.  https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/bam And no, samtools fastq will not make paired end fastqs suitable to rerun through cellranger, because the read data from the read containing the cell barcode and UMI will not have its own line in the bams.

@swbarnes2 There are indeed many tags, but isn't the UMI tag different than UMI counts? So how would I derive the expression matrix from these said tags? Thank you btw for your answer!

Each read has its UMI as a bam tag.  You can go through the bam line by line and figure out which reads were assigned to what cell barcode, UMI and gene.

This reverse process is also complicated by the fact that error-correction is performed on the barcodes and UMIs as well. For scRNA-seq, you typically start with 4 FASTQs: UMIs, barcodes, R1, and R2. There will likely be errors in those sequences, so barcodes/UMIs are corrected for by comparing to a whitelist of known input barcodes. You can take a look at the "Feature Barcoding Tags" section of the 10x support page @swbarnes2 listed to see what tags to look for. As mentioned, it's not a simple BAM -> FASTQ operation, but you should be able to do it

@James, that's not at all how the fastqs come out for 10xGenomics samples. (Or dropSeq either)  And feature barcoding is something new and extra; an old experiment probably won't have it.

My mistake, I was mixing up the scATAC-seq and scRNA-seq protocols. The scATAC protocol will have 3 or 4 FASTQs (https://support.10xgenomics.com/single-cell-atac/software/pipelines/latest/using/fastq-input), whereas scRNA has 2 or 3 (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/fastq-input). They both do barcode correction, but the input FASTQs are slightly different, and you're correct that the UMIs are probably only on newer datasets

@James UMIs have been used from at least cellranger 1.1.0 in 2016, if not earlier. Features are different, and were just introduced this year.

Cell Ranger bam files have the corrected and uncorrected barcodes and UMIs.  It is definitely possible to go from a Cell Ranger bam file to fastqs. 10X has a tool for doing it (https://support.10xgenomics.com/docs/bamtofastq)


text:
expansions and contractions from OrthoFinder
The input file looks like this, and the complete file can be found here:
OG0000008 NbL07g11380.1 NbL19g07810.1 NbL19g09170.1 NbL19g19070.1 NbQ01g01670.1 NbQ01g03330.1 NbQ01g04070.1 NbQ01g04670.1 NbQ01g05120.1 NbQ01g05870.1 NbQ01g06940.1 NbQ01g07580.1 NbQ01g08860.1 NbQ01g10050.1 NbQ01g10360.1 NbQ01g14200.1 NbQ01g14790.1 NbQ01g16080.1 NbQ01g17760.1 NbQ01g19270.1 NbQ01g19310.1 NbQ01g19390.1 NbQ01g21260.1 NbQ01g21330.1 NbQ01g21740.1 NbQ01g21910.1 NbQ01g23100.1 NbQ01g24620.1 NbQ01g25340.1 NbQ01g26060.1 NbQ01g26320.1 NbQ02g00750.1 NbQ02g03100.1 NbQ02g03420.1 NbQ02g03610.1 NbQ02g03680.1 NbQ02g05120.1 NbQ02g07460.1 NbQ02g08170.1 NbQ02g08330.1 NbQ02g09220.1 NbQ02g09400.1 NbQ02g10620.1 NbQ02g11310.1 NbQ02g14330.1 NbQ02g14460.1 NbQ02g14520.1 NbQ02g15320.1 NbQ02g17090.1 NbQ02g17130.1 NbQ02g20290.1 NbQ02g23070.1 NbQ02g23420.1 NbQ02g24450.1 NbQ02g24480.1 NbQ02g26700.1 NbQ03g00830.1 NbQ03g01970.1 NbQ03g04460.1 NbQ03g06900.1 NbQ03g09530.1 NbQ03g10620.1 NbQ03g12760.1 NbQ03g13450.1 NbQ03g15540.1 NbQ03g15640.1 NbQ03g17180.1 NbQ03g20740.1 NbQ03g21510.1 NbQ03g24670.1 NbQ04g01350.1 NbQ04g01720.1 NbQ04g08420.1 NbQ04g09090.1 NbQ04g10450.1 NbQ04g11470.1 NbQ04g12120.1 NbQ04g14130.1 NbQ04g15440.1 NbQ04g15860.1 NbQ04g16450.1 NbQ04g16620.1 NbQ04g17760.1 NbQ04g19040.1 NbQ04g20020.1 NbQ05g03320.1 NbQ05g04660.1 NbQ05g05970.1 NbQ05g07500.1 NbQ05g08900.1 NbQ05g09760.1 NbQ05g10830.1 NbQ05g11150.1 NbQ05g11340.1 NbQ05g11510.1 NbQ05g11530.1 NbQ05g11780.1 NbQ05g16980.1 NbQ05g18190.1 NbQ05g21710.1 NbQ05g23400.1 NbQ06g01110.1 NbQ06g01430.1 NbQ06g04200.1 NbQ06g04440.1 NbQ06g05330.1 NbQ06g05770.1 NbQ06g05820.1 NbQ06g06700.1 NbQ06g08620.1 NbQ06g09190.1 NbQ06g10460.1 NbQ06g15220.1 NbQ06g15330.1 NbQ06g15700.1 NbQ06g16320.1 NbQ06g16590.1 NbQ06g17590.1 NbQ06g17670.1 NbQ06g20050.1 NbQ07g01030.1 NbQ07g02010.1 NbQ07g04350.1 NbQ07g04900.1 NbQ07g05610.1 NbQ07g06200.1 NbQ07g07110.1 NbQ07g07690.1 NbQ07g08640.1 NbQ07g10390.1 NbQ07g11920.1 NbQ07g14130.1 NbQ07g15590.1 NbQ07g15620.1 NbQ07g16910.1 NbQ07g17130.1 NbQ07g17950.1 NbQ08g00060.1 NbQ08g02240.1 NbQ08g02300.1 NbQ08g02310.1 NbQ08g03290.1 NbQ08g05330.1 NbQ08g09280.1 NbQ08g14890.1 NbQ08g15820.1 NbQ08g15950.1 NbQ08g19830.1 NbQ08g20150.1 NbQ08g22050.1 NbQ08g22620.1 NbQ09g02100.1 NbQ09g02620.1 NbQ09g03950.1 NbQ09g04200.1 NbQ09g06040.1 NbQ09g06640.1 NbQ09g08160.1 NbQ09g08330.1 NbQ09g09660.1 NbQ09g11220.1 NbQ09g13860.1 NbQ09g15180.1 NbQ09g15310.1 NbQ09g16530.1 NbQ09g17900.1 NbQ09g18100.1 NbQ09g18720.1 NbQ09g19280.1 NbQ09g21840.1 NbQ10g00480.1 NbQ10g01350.1 NbQ10g02870.1 NbQ10g03640.1 NbQ10g03730.1 NbQ10g08070.1 NbQ10g09510.1 NbQ10g11010.1 NbQ10g11760.1 NbQ10g12050.1 NbQ10g12060.1 NbQ10g12910.1 NbQ10g19200.1 NbQ10g19930.1 NbQ10g20390.1 NbQ10g20730.1 NbQ10g21080.1 NbQ10g21140.1 NbQ10g24010.1 NbQ11g00310.1 NbQ11g01210.1 NbQ11g01370.1 NbQ11g04610.1 NbQ11g04800.1 NbQ11g06060.1 NbQ11g07820.1 NbQ11g08390.1 NbQ11g09100.1 NbQ11g09350.1 NbQ11g13660.1 NbQ11g13930.1 NbQ11g16260.1 NbQ11g17360.1 NbQ11g18430.1 NbQ11g21080.1 NbQ11g23280.1 NbQ11g23990.1 NbQ11g25050.1 NbQ12g03770.1 NbQ12g04850.1 NbQ12g07340.1 NbQ12g09080.1 NbQ12g10820.1 NbQ12g12070.1 NbQ12g14750.1 NbQ12g15000.1 NbQ12g15230.1 NbQ12g20380.1 NbQ12g21080.1 NbQ12g21830.1 NbQ12g23960.1 NbQ13g01300.1 NbQ13g02350.1 NbQ13g03860.1 NbQ13g04410.1 NbQ13g08800.1 NbQ13g09850.1 NbQ13g10370.1 NbQ13g11700.1 NbQ13g12420.1 NbQ13g15780.1 NbQ13g16040.1 NbQ13g23160.1 NbQ13g24120.1 NbQ13g24540.1 NbQ13g25080.1 NbQ13g25490.1 NbQ13g28240.1 NbQ13g29770.1 NbQ14g01070.1 NbQ14g03950.1 NbQ14g05360.1 NbQ14g05410.1 NbQ14g06880.1 NbQ14g07270.1 NbQ14g07500.1 NbQ14g10290.1 NbQ14g10770.1 NbQ14g14320.1 NbQ14g17890.1 NbQ14g18710.1 NbQ14g20960.1 NbQ14g22890.1 NbQ15g00150.1 NbQ15g02300.1 NbQ15g02330.1 NbQ15g02350.1 NbQ15g03230.1 NbQ15g06190.1 NbQ15g07120.1 NbQ15g07750.1 NbQ15g09000.1 NbQ15g09050.1 NbQ15g11920.1 NbQ15g12650.1 NbQ15g12840.1 NbQ15g15670.1 NbQ15g15930.1 NbQ15g18670.1 NbQ15g19070.1 NbQ15g20620.1 NbQ15g22880.1 NbQ15g23000.1 NbQ15g26060.1 NbQ16g00880.1 NbQ16g04360.1 NbQ16g06490.1 NbQ16g09100.1 NbQ16g11020.1 NbQ16g11560.1 NbQ16g13810.1 NbQ16g13820.1 NbQ16g17040.1 NbQ16g17130.1 NbQ16g17340.1 NbQ16g18390.1 NbQ16g18430.1 NbQ16g23100.1 NbQ16g23570.1 NbQ16g24270.1 NbQ16g25200.1 NbQ16g25830.1 NbQ16g25880.1 NbQ16g25990.1 NbQ16g26610.1 NbQ16g26660.1 NbQ16g28010.1 NbQ16g28180.1 NbQ17g01150.1 NbQ17g01180.1 NbQ17g01570.1 NbQ17g01950.1 NbQ17g05460.1 NbQ17g05540.1 NbQ17g05980.1 NbQ17g07990.1 NbQ17g08300.1 NbQ17g09330.1 NbQ17g09400.1 NbQ17g10090.1 NbQ17g11220.1 NbQ17g13030.1 NbQ17g15460.1 NbQ17g16690.1 NbQ17g20980.1 NbQ17g22370.1 NbQ17g25040.1 NbQ17g28730.1 NbQ18g02140.1 NbQ18g02740.1 NbQ18g05440.1 NbQ18g06120.1 NbQ18g07470.1 NbQ18g12320.1 NbQ18g12530.1 NbQ18g12850.1 NbQ18g13840.1 NbQ18g14420.1 NbQ18g14930.1 NbQ18g15730.1 NbQ18g17750.1 NbQ18g17850.1 NbQ18g21060.1 NbQ19g01040.1 NbQ19g05480.1 NbQ19g06450.1 NbQ19g06510.1 NbQ19g08330.1 NbQ19g11840.1 NbQ19g11880.1 NbQ19g13750.1 NbQ19g14190.1 NbQ19g14210.1 NbQ19g14920.1 NbQ19g18540.1 NbQ19g19870.1 NbQ19g21020.1 NbQ19g21220.1 NbQ19g22080.1 NbQ19g22800.1 NbQ19g24690.1 NbQ19g24730.1 rna19561
OG0000001 Capann_59V1aChr01g048170.1 NbL01g00020.1 NbL01g00940.1 NbL01g02330.1 NbL01g03550.1 NbL01g03650.1 NbL01g04410.1 NbL01g04920.1 NbL01g06850.1 NbL01g16120.1 NbL01g19150.1 NbL01g20140.1 NbL01g20930.1 NbL01g22230.1 NbL01g24190.1 NbL01g24280.1 NbL01g24300.1 NbL02g00570.1 NbL02g00900.1 NbL02g01270.1 NbL02g02110.1 NbL02g02210.1 NbL02g02470.1 NbL02g03180.1 NbL02g04740.1 NbL02g04750.1 NbL02g06120.1 NbL02g06860.1 NbL02g07280.1 NbL02g07680.1 NbL02g07740.1 NbL02g09780.1 NbL02g11320.1 NbL02g12670.1 NbL02g13080.1 NbL02g14050.1 NbL02g14190.1 NbL02g15010.1 NbL02g15890.1 NbL02g16190.1 NbL02g16730.1 NbL02g17070.1 NbL02g17360.1 NbL02g18820.1 NbL02g19340.1 NbL02g20100.1 NbL02g23950.1 NbL02g24800.1 NbL03g01610.1 NbL03g01680.1 NbL03g01890.1 NbL03g02230.1 NbL03g02600.1 NbL03g03410.1 NbL03g04990.1 NbL03g05400.1 NbL03g08030.1 NbL03g08250.1 NbL03g08690.1 NbL03g10230.1 NbL03g11060.1 NbL03g13030.1 NbL03g14960.1 NbL03g15110.1 NbL03g16690.1 NbL03g16900.1 NbL03g18260.1 NbL03g18950.1 NbL03g21180.1 NbL03g21210.1 NbL03g21530.1 NbL03g22960.1 NbL03g24430.1 NbL04g01140.1 NbL04g01490.1 NbL04g02030.1 NbL04g02560.1 NbL04g03700.1 NbL04g04160.1 NbL04g05240.1 NbL04g05420.1 NbL04g05850.1 NbL04g12420.1 NbL04g12640.1 NbL04g13650.1 NbL04g13780.1 NbL04g14310.1 NbL04g16260.1 NbL04g17750.1 NbL04g18380.1 NbL04g18870.1 NbL04g19030.1 NbL04g19630.1 NbL05g00320.1 NbL05g03060.1 NbL05g03300.1 NbL05g04060.1 NbL05g07620.1 NbL05g08630.1 NbL05g09580.1 NbL05g10060.1 NbL05g11400.1 NbL05g12280.1 NbL05g13170.1 NbL05g16020.1 NbL05g17530.1 NbL05g17730.1 NbL05g18340.1 NbL05g18590.1 NbL05g18600.1 NbL05g19640.1 NbL05g20640.1 NbL05g21000.1 NbL05g21640.1 NbL05g22610.1 NbL06g00640.1 NbL06g00660.1 NbL06g02210.1 NbL06g03150.1 NbL06g03680.1 NbL06g04910.1 NbL06g07950.1 NbL06g09970.1 NbL06g11480.1 NbL06g12220.1 NbL06g12400.1 NbL06g12460.1 NbL06g12850.1 NbL06g13120.1 NbL06g14450.1 NbL06g14780.1 NbL06g16990.1 NbL06g17200.1 NbL06g17760.1 NbL06g20380.1 NbL07g02100.1 NbL07g02540.1 NbL07g02970.1 NbL07g03110.1 NbL07g04840.1 NbL07g05350.1 NbL07g06580.1 NbL07g07530.1 NbL07g08450.1 NbL07g09380.1 NbL07g09870.1 NbL07g10730.1 NbL07g10850.1 NbL07g11080.1 NbL07g12450.1 NbL07g12710.1 NbL07g13110.1 NbL07g13920.1 NbL07g14240.1 NbL07g15520.1 NbL07g16220.1 NbL07g17480.1 NbL08g01820.1 NbL08g02750.1 NbL08g02930.1 NbL08g03510.1 NbL08g03620.1 NbL08g03850.1 NbL08g03970.1 NbL08g04040.1 NbL08g06150.1 NbL08g06410.1 NbL08g06680.1 NbL08g06730.1 NbL08g07620.1 NbL08g08450.1 NbL08g08640.1 NbL08g09910.1 NbL08g10160.1 NbL08g11760.1 NbL08g12570.1 NbL08g13630.1 NbL08g13890.1 NbL08g15050.1 NbL08g15340.1 NbL08g18010.1 NbL08g18420.1 NbL08g19080.1 NbL08g19190.1 NbL09g00900.1 NbL09g02160.1 NbL09g02330.1 NbL09g02470.1 NbL09g04150.1 NbL09g05210.1 NbL09g07010.1 NbL09g09070.1 NbL09g10290.1 NbL09g10500.1 NbL09g11220.1 NbL09g13490.1 NbL09g15290.1 NbL09g15830.1 NbL09g17240.1 NbL09g19250.1 NbL09g19460.1 NbL09g20190.1 NbL09g21040.1 NbL09g21520.1 NbL09g23460.1 NbL10g00080.1 NbL10g03710.1 NbL10g04330.1 NbL10g04560.1 NbL10g05200.1 NbL10g06320.1 NbL10g07510.1 NbL10g07960.1 NbL10g08670.1 NbL10g08970.1 NbL10g11120.1 NbL10g11340.1 NbL10g11820.1 NbL10g13720.1 NbL10g14560.1 NbL10g14770.1 NbL10g16430.1 NbL10g18140.1 NbL10g18380.1 NbL10g19280.1 NbL10g19690.1 NbL10g21210.1 NbL10g22680.1 NbL10g23160.1 NbL10g23560.1 NbL10g24210.1 NbL11g00680.1 NbL11g00970.1 NbL11g01230.1 NbL11g01270.1 NbL11g01520.1 NbL11g01530.1 NbL11g02920.1 NbL11g03540.1 NbL11g03990.1 NbL11g05630.1 NbL11g08950.1 NbL11g08980.1 NbL11g09510.1 NbL11g10840.1 NbL11g11030.1 NbL11g11230.1 NbL11g12430.1 NbL11g13300.1 NbL11g15430.1 NbL11g16390.1 NbL11g16410.1 NbL11g17320.1 NbL11g18090.1 NbL11g21310.1 NbL11g21470.1 NbL11g21780.1 NbL11g21820.1 NbL11g22270.1 NbL11g22310.1 NbL11g23180.1 NbL11g24100.1 NbL12g00130.1 NbL12g01810.1 NbL12g02230.1 NbL12g02720.1 NbL12g02760.1 NbL12g04120.1 NbL12g04550.1 NbL12g06630.1 NbL12g07830.1 NbL12g09170.1 NbL12g10580.1 NbL12g12090.1 NbL12g12490.1 NbL12g12630.1 NbL12g12800.1 NbL12g13320.1 NbL12g13460.1 NbL12g14430.1 NbL12g14970.1 NbL12g15490.1 NbL12g17460.1 NbL12g18190.1 NbL12g18590.1 NbL12g19900.1 NbL12g20690.1 NbL12g22040.1 NbL12g22560.1 NbL13g00350.1 NbL13g01440.1 NbL13g02400.1 NbL13g03210.1 NbL13g03360.1 NbL13g04070.1 NbL13g05250.1 NbL13g08460.1 NbL13g09010.1 NbL13g09140.1 NbL13g10290.1 NbL13g11570.1 NbL13g13370.1 NbL13g14910.1 NbL13g18680.1 NbL13g19510.1 NbL13g23520.1 NbL13g24010.1 NbL13g24190.1 NbL13g24460.1 NbL13g26310.1 NbL13g26640.1 NbL13g26860.1 NbL13g27260.1 NbL13g27960.1 NbL14g02460.1 NbL14g02750.1 NbL14g08750.1 NbL14g08910.1 NbL14g09120.1 NbL14g09540.1 NbL14g09920.1 NbL14g11070.1 NbL14g11150.1 NbL14g12570.1 NbL14g14530.1 NbL14g14860.1 NbL14g15240.1 NbL14g15460.1 NbL14g16620.1 NbL14g16910.1 NbL14g17940.1 NbL14g21150.1 NbL14g21750.1 NbL14g21910.1 NbL15g00790.1 NbL15g01170.1 NbL15g02310.1 NbL15g04220.1 NbL15g05970.1 NbL15g06340.1 NbL15g06440.1 NbL15g07020.1 NbL15g07370.1 NbL15g07470.1 NbL15g09010.1 NbL15g13210.1 NbL15g14550.1 NbL15g14600.1 NbL15g17290.1 NbL15g18170.1 NbL15g19710.1 NbL15g21840.1 NbL15g21930.1 NbL15g23410.1 NbL15g23420.1 NbL15g23430.1 NbL15g25130.1 NbL15g25200.1 NbL16g01760.1 NbL16g02140.1 NbL16g04460.1 NbL16g05010.1 NbL16g05020.1 NbL16g06780.1 NbL16g07540.1 NbL16g07980.1 NbL16g09760.1 NbL16g10610.1 NbL16g12320.1 NbL16g13510.1 NbL16g14420.1 NbL16g15690.1 NbL16g17420.1 NbL16g17790.1 NbL16g17880.1 NbL16g18730.1 NbL16g18940.1 NbL16g19440.1 NbL16g20980.1 NbL16g23180.1 NbL16g23610.1 NbL16g23660.1 NbL16g23910.1 NbL16g24550.1 NbL16g24640.1 NbL16g25300.1 NbL16g25630.1 NbL16g26710.1 NbL17g01590.1 NbL17g02070.1 NbL17g02120.1 NbL17g02920.1 NbL17g03040.1 NbL17g03540.1 NbL17g03700.1 NbL17g03800.1 NbL17g05400.1 NbL17g07510.1 NbL17g08450.1 NbL17g08930.1 NbL17g10090.1 NbL17g14370.1 NbL17g14600.1 NbL17g15390.1 NbL17g15900.1 NbL17g16000.1 NbL17g16910.1 NbL17g17480.1 NbL17g18240.1 NbL17g20020.1 NbL17g20830.1 NbL17g21220.1 NbL17g21690.1 NbL17g25960.1 NbL18g00030.1 NbL18g00150.1 NbL18g00310.1 NbL18g00670.1 NbL18g00700.1 NbL18g01630.1 NbL18g02650.1 NbL18g04460.1 NbL18g05210.1 NbL18g05690.1 NbL18g07270.1 NbL18g07440.1 NbL18g07500.1 NbL18g09090.1 NbL18g09810.1 NbL18g09880.1 NbL18g10500.1 NbL18g10990.1 NbL18g12070.1 NbL18g13060.1 NbL18g17480.1 NbL19g00230.1 NbL19g04470.1 NbL19g04700.1 NbL19g07770.1 NbL19g07870.1 NbL19g09260.1 NbL19g09870.1 NbL19g13480.1 NbL19g14360.1 NbL19g14400.1 NbL19g14720.1 NbL19g17700.1 NbL19g18860.1 NbL19g22750.1 NbQ13g00370.1 NbQ14g10310.1 NbQ17g06050.1

This script appears to have some problems with the file:
import pandas as pd

orthofinder_output = "../OrthoFinder-res/Results_Jan23/Orthogroups/Orthogroups-fixed.txt"
orthogroups = pd.read_csv(orthofinder_output, sep=' ', header=None)

# Extract gene family information
expansions = {}
contractions = {}
for i, row in orthogroups.iterrows():
    if len(row) > 2:
        expansions[row[0]] = row[2:]
    else:
        contractions[row[0]] = row[1]

# Print expansions and contractions
print("Expansions:", expansions)
print("Contractions:", contractions)

I got the following error:
geneFamilyExpansionsContractions.py:
 20: DtypeWarning:
   Columns (2,3,4,5,6,...,964,965,966,967,968) have mixed types.
   Specify dtype option on import or set low_memory=False.
  orthogroups = pd.read_csv(orthofinder_output, sep=' ', header=None)
...
Expansions:
...
Name: 63241, Length: 967, dtype: object, 'OG0063242':
2      NaN
3      NaN
4      NaN
5      NaN
6      NaN
      ...
964    NaN
965    NaN
966    NaN
967    NaN
968    NaN
Name: 63242, Length: 967, dtype: object}
Contractions: {}

How can it be fixed?

The answer is actually right there in the error message (warning message, actually):

DtypeWarning: Columns (2,... 968) have mixed types. Specify dtype option on import or set low_memory=False.

Setting low_memory-False makes the error go away:
orthogroups = pd.read_csv(orthofinder_output, sep=' ', header=None, low_memory=False)

According to this random blog post I found when I searched for "pandas DtypeWarning: have mixed types", this is because:

Pandas is really nice because instead of stopping altogether, it guesses which dtype a column has. However, by default, Pandas has the low_memory=True argument. This means that the CSV file gets split up in multiple chunks and the guess is done for every chunk, resulting in a column with multiple dtypes. Pandas is so kind to let you know it was confused and something might have happened.

The same blog also suggests that this isn't really the right approach, and that makes sense, since this is forcing Pandas to read the entire file and that might be a problem for huge files or low memory machines. So, reading on, the blog suggests the other option provided by the error, setting the dtype manually.
Now, I am afraid I have never used Pandas and don't really know much Python either, but a quick glance at the Pandas documentation suggests that the object data type is the most versatile and should work for everything, including strings. Since your file is nothing but a collection of strings, as far as I can tell, this seems like the best approach.
So I changed your script to this:
#!/usr/bin/env python3

import pandas as pd
import sys

orthofinder_output = sys.argv[1]
orthogroups = pd.read_csv(orthofinder_output, sep=' ', header=None, low_memory=False)
# Extract gene family information
expansions = {}
contractions = {}
for i, row in orthogroups.iterrows():
    if len(row) > 2:
        expansions[row[0]] = row[2:]
    else:
        contractions[row[0]] = row[1]

# Print expansions and contractions
print("Expansions:", expansions)
print("Contractions:", contractions)

And it seems to work just fine with no warnings shown.

The error right at the end is

Specify dtype option on import or set low_memory=False.

Edit Oh wait @terdon has already spotted this! :-D ... great minds and all that :)
In addition, the pandas print out it states ...
dtype object

This isn't good. Name: 63242, is a bit unusual BTW worth checking.
Looking at the values ...
NaN

... that is really not good.
What I think is happening is it's been assigned to a numerical datatype, but it's categorical data. Hence NaN because like the alphabet can be none-base 10 numbers, in data compression used a lot, but it's got to be consistent pattern (I genuinely don't know if pandas can spot/work in none-base 10 number systems). pandas is really C (probably C++) - its a static type backend - so those datatypes are the key to this.
It makes sense because your data is categorical data and NaN means its managed to confuse a large amount of it for e.g. int64.
Looking the documentation the alternative is
orthogroups = pd.read_csv(orthofinder_output, sep=' ', header=None, dtype='category')

The documentation is

low_memorybool, default True
Internally process the file in chunks, resulting in lower memory use while parsing, but possibly mixed type inference. To ensure no mixed types either set False, or specify the type with the dtype parameter. Note that the entire file is read into a single DataFrame regardless, use the chunksize or iterator parameter to return the data in chunks. (Only valid with C parser).

What you want is to issue a df.dtypes command once it is loaded and the output will define the datatype per column,
'A'    int64
'B'    int32
'C'   category

If you have mixed dtypes it may (I stress may) be able to correct it via,
df.astype('int64') # this is not your data
df.astype('category') # this is your data

HOWEVER, I don't think this will work, because NaN indicates its trashed the data. Anyway the data is definitely categorical BTW.
Anyway, good question and interesting stuff.

As a bit of gossip, you can also switch to categorical ordered data I think this will work for dataframes, this example is Series. There's other ways to order the data BTW.
from pandas.api.types import CategoricalDtype
my_dtype = CategoricalDtype(\
    categories=[2, 1], ordered=True)
df.dtype(my_dtype)

The output is 1 and 2 now switch places because 2 < 1. Categorial datatypes look to be the right datatype for your inputs.


text:
SNP Signatures with Limited WGS Data:
On 80 WGS samples, I'm dissecting SNP signatures linked to milk production in a scarcely studied animal. Post-variant calling and QC association analysis have been tricky. I'm here to tap into our collective wisdom and refine my approach.
Study Snapshot:
Objective: Identify SNPs that differentiate between high (>5 L/day) and low (<5 L/day) milk producers.
Data: 80 WGS samples.
Methodology:
Phenotypes were normalized (high = 1, low = 0).
PCs from PCA are used to correct for population stratification, with phenotype information as a covariate.
Initial approach: linear models, then FDR and Bonferroni for multiple comparison correction. A significant SNP threshold was set at <1e-5, yielding four significant SNPs.
Expanded to Linear Mixed Models (LMM), we found six significant SNPs at 1e-5 and 99 at 1e-4, focusing on positive associations with high milk production.
Seeking Guidance:
Feasibility: Can 80 WGS samples provide a reliable genetic signature for milk production?
Methodology Recommendations: Given my data and initial results, what methodologies or modifications could improve the analysis?

Visualise your problem and solution.
6 SNPs is not a particularly large search space (64 possible combinations assuming a dominant or recessive model with dimorphisms, 729 if heterozygous variants are taken into account). I would graph all those haplotypes out in a heatmap (but only the ones present in your samples, so there'd be a maximum of 80 different ones) and make sure that there is a clear pattern.
If population stratification is a concern, then you should also show the PCA, with separate colouring for the different phenotypes. If the two groups have similar spread in the PCA (which would be expected for a trait that is only affected by a small number of genes), then PCA correction doesn't make sense (and may lead to false associations).


text:
How to create a legend in FeaturePlot with do.hover at the same time?
you can see in my feature plot I can only see the cells from which sample but I cannot see the legend. And if I change the command I can see the legend but cannot see the sample identities.FeaturePlot(object = VarMerge, features.plot = c("Hcn4"), do.hover = T, no.legend = F,)
if I change the command to
FeaturePlot(object = VarMerge, features.plot = c("Hcn4"), no.legend = F

it will show like the second picture. So I just want to know whether in the first picture can I do something to the command and add the legend like in second picture?


text:
Parsing gene expression value for a list of genes from GEO
I am looking for the gene expression data in GEO GSE15543 with the following keywords,
(Pancreatic Beta Cells) AND "Homo sapiens"[porgn:__txid9606]

I would like to parse the gene expression value of a list of genes from all the samples reported in any particular study.I understand there are two different formats for downloading datasets, SOFT and MINiML.
I have downloaded the soft.gz file and I am using BioConductor to read the contents from soft file using the instructions given here . The file has the content that I 'm looking for , I have checked the contents on the geo2r webpage. I'm still facing challenge in finding out how to fetch out the gene expression values for a list of genes. Any suggestions?
From the files that I download, I wish to parse the gene expression value and the Sample Title and Sample Description so that I can classify the values based on disease state and control state. Any suggestions on which format can be downloaded for the above-mentioned task?

GEO provides a GUI for this (just click on the "Analyze with GEO2 button"), you can even check/download the R/Bioconductor code used for the analyses.
Basically you will create groups of samples to compare/contrast and R/Bioconductor will run behind the scenes. Here is a nice tutorial.


text:
How to detail the specific GO terms
How can I get more specific GO terms when using  clusterProfiler? I got my Dotplot by:
cluster0_all_gene_ego<-enrichGO(cluster0_all_gene,keyType = 'SYMBOL',
                                OrgDb = org.Mm.eg.db,
                                ont = "BP",
                                pAdjustMethod = "BH",
                                pvalueCutoff = 1,
                                qvalueCutoff = 1)
dotplot(cluster0_all_gene_ego)

as below:
It seems that I can only get the parent terms as some one tell me that there can have more specific child terms like below:
GO:0008150 biological_process
     ￼ GO:0048856 anatomical structure development
     ￼ GO:0032501 multicellular organismal process
         ￼ GO:0032502 developmental process
         ￼ GO:0007275 multicellular organism development
             ￼ GO:0048731 system development
             ￼ GO:0009888 tissue development
                 ￼ GO:0048513 animal organ development
                 ￼ GO:0072359 circulatory system development
                 ￼ GO:0060537 muscle tissue development
                     ￼ GO:0048738 cardiac muscle tissue development
                     ￼ GO:0007507 heart development
                     ￼ GO:0014706 striated muscle tissue development
                         ￼ GO:0003161 cardiac conduction system development
                             ￼ GO:0003163 sinoatrial node development
                                 ￼ GO:0060921 sinoatrial node cell differentiation

I would like to have some terms like the sinoatrial node, so did anyone know how to modify the code?
I have got the GO term's result and find most of the gene is expressed in the sinoatrial node, but there are no sinoatrial node terms.


text:
Understanding fast5 format
I would like to reproduce the results on figure 1 from THIS publication.
I downloaded the data and explore it with h5ls and h5dump based on HDF5.
I would like to understand:

h5ls -r <_fast5file_> | grep Dataset:

/Analyses/Basecall_2D_000/BaseCalled_2D/Alignment Dataset {6925}
/Analyses/Basecall_2D_000/BaseCalled_2D/Fastq Dataset {SCALAR}
/Analyses/Basecall_2D_000/BaseCalled_complement/Events Dataset {5576}
/Analyses/Basecall_2D_000/BaseCalled_complement/Fastq Dataset {SCALAR}
/Analyses/Basecall_2D_000/BaseCalled_complement/Model Dataset {4096}
/Analyses/Basecall_2D_000/BaseCalled_template/Events Dataset {5793}
/Analyses/Basecall_2D_000/BaseCalled_template/Fastq Dataset {SCALAR}
/Analyses/Basecall_2D_000/BaseCalled_template/Model Dataset {4096}
/Analyses/Basecall_2D_000/HairpinAlign/Alignment Dataset {5107}
/Analyses/Basecall_2D_000/Log Dataset {SCALAR}
/Analyses/Calibration_Strand_000/Log Dataset {SCALAR}
/Analyses/EventDetection_000/Reads/Read_0/Events Dataset {11465}

Alignment? Then complement/Fastq is not the complement of template/Fastq (analyzed with h5dump). What is it? Also, events have sometimes 4

Some other datasets that were acquired somewhere else have the raw data and fastq but there is no associated information about kmer and mean(current)/std/length/etc... How do I generate events from raw data?

I went through all documentations and tutorials that I could found in the web but couldn't find anything simple and complete. It feels as if it was a secret.
Example fast5file would be ecoli_er2925.MSssI.timp.100215.fast5/pass/kelvin_100215_ecolisssI_2515_1_ch195_file0_strand.fast5 from PRJEB13021

FAST5 is a proprietary format developed by Oxford Nanopore Technologies, and there's not much great documentation online. All FAST5 files will have the Raw/ field, which contains the original measured Signal. Additional Analyses/ fields can be added by tools such as basecallers which convert signal->ACGT (Guppy) or signal-level aligners (Tombo).
In particular, /Analyses/Basecall_2D_000/ contains the results of 2D basecalling, added by Guppy when run with the --post-out flag. I think Tombo adds the event-level annotations you mention (k-mer mean/std/length). In 2D basecalling, both the physical forward and reverse strands for the same stretch of DNA are measured and basecalled. The basecalls will contain errors, so these will not be exact reverse complements of one another.
Hope this answers your questions, I have a slightly more in-depth explanation on my blog: https://timd.one/blog/genomics/file_formats.php

FWIW, this particular FAST5 was generated in 2017, before Guppy was available

I'll give a different perspective. Oxford Nanopore FAST5 is NGS within a HDF5 database. A clear description is here.
HDF5 is a format used for rapid I/O of large data sets.
Oxford Nano appear to have written a bit boilerplate around HDF5 and called it something that sounds like 'fasta'.
A cool part of HDF5 (not FAST5) is  HDF5 View, which can be downloaded here. Its cool because you can visualise the file structure of HDF5 db via a GUI. I very much suspect if you use HDF5 View on your FAST5 database you will be able to navigate (see not change) the entire file structure.
To answer your question directly h5ls is listing the HDF5 database. To construct a HDF5 db you assign the data under a series of 'keys' which forms the directory structure. However, you perhaps might not need Oxford tools, because  Python pandas have direct series of methods to connect, construct and manipulate HDF5.
I raised the HDF5 in a question then answered here on Bioinfo SE here.
Again the reason for HDF5's importance is unlike an SQL connection, data transfer and disconnect, HDF5 shifts data very quickly.
I do question how Oxford Nanopore can lay claim to FAST5 because HDF5 is open source and always has been. If they want to lay claim to some boilerplate, be my guest - I write the stuff all the time ;-)
HDF5 however is a serious project with major development time and appears to still be under active development.


text:
How to load output from the kegg API in biopython into a pandas dataframe?
Biopython provides a (bit unintuitive) API to access to the Kyoto Encyclopedia of Genes and Genomes (KEGG). I am trying to make use of it, but the output is quite unhandy as a string. What is the best way to parse the data into a pandas.DataFrame()?
from Bio.KEGG import REST as Kegg
import io

pd.read_table(kegg.kegg_get('K02545'))

solution from @sören:
I came up with this solution for my problem:
from Bio.KEGG import REST as kegg

def _get_kegg(kegg_id):
    kegg_output = kegg.kegg_get(kegg_id).read()
    results = {}
    for line in kegg_output.split('\n'):
        splits = line.split()
        if not line.startswith(' '):
            if len(splits) > 0:
                key = splits[0]
                value = ' '.join(splits[1:])
                results[key] = value
        else:
            results[key] += ' '.join(splits)
    return pd.DataFrame(results, index=[kegg_id])

_get_kegg_v = np.vectorize(_get_kegg)

def get_kegg_info(kegg_ids):
    if isinstance(kegg_ids, str):
        kegg_ids = [kegg_ids]
    return pd.concat(_get_kegg_v(kegg_ids), sort=False)


text:
RNA seq fasta file annotation from alignment to reference matches
I've got a fasta file with some RNA seq data and another csv file with the output from plast where I've aligned it to a reference using plastn. I'm struggling with figuring out a command to append my original RNA seq data fasta headers with the annotations.
I've looked at awk and sed but can't wrap my head around how to use them in this instance.
CSV output example:
Contig11    ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA   78.95   95  20  0   371 465 4627    4721    1e-14   82.4    75  625 1   0   15.2    0   5692    1   0   15.2    0

And example fasta contig:
>Contig11
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA

How would I be able to parse through the csv to find a matching contig with the fasta file, and append the second csv column to the fasta header. Ideally to produce:
>Contig11 ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA....

Cheers for any help!

What format is your "csv" in? It isn't a "real" csv (that means "comma separated file" where columns are defined by commas), so what is the column separator? Is it a space or a tab?

An awk approach (assuming your input "csv" file is space-separated):
$ awk '(NR==FNR){a[">"$1]=$2; next;} (/>/){$1=$1" "a[$1]}1;' foo.csv foo.fa
>Contig11 ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA

Explanation

NR==FNR : FNR is the line number of the current input file and NR is the line number of all input. When processing more than one file, the two will only be equal while reading the first file.
(NR==FNR){a[">"$1]=$2; next;} : if we're reading the first file, the csv, save the second field as the value in the associative array a, whose indices are a > followed by the 1st field. I use the > only to simplify matching it later in the fasta file. The next moves on to the next line so the rest of the command isn't run for the 1st file.
(/>/){$1=$1" "a[$1]}: Now we're in the 2nd file. If the line starts with a >, then set the 1st field to be itself, a space and then whatever was stored in the array a for this 1st field.
1; : this is awk shorthand for "print this line".

You could expand the above into:
awk '{
        if(NR==FNR){
            a[">"$1]=$2;
        }
        else{
            if(/>/){
                $1=$1" "a[$1]
            }
        }
        print;
      }' foo.csv foo.fa

If your files are small enough to be sorted easily, you can use the FastaToTbl and TblToFasta scripts I have posted before and then join:
$ join -o 1.1,2.2,1.2 <(FastaToTbl foo.fa | sort ) <(sort  foo.csv ) | TblToFasta
>Contig11 ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA
>Contig12 mamamama|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
AKHSALKJHDLKAJBDLKAJDBADADKJALDKJBAD


text:
convert supplementary reads to primary in sam or bam
I have a problem with a tool, that possibly ignores supplementary reads. I want to find out a bit how this tool works by converting all supplementary reads to primary reads, and then change the names of those reads (because I guess that primary reads with the same name are ignored?). If the tool produces something, instead of nothing (like now), then I know for sure its the supplementary read flag that could be the problem.
I use SV detection tools. The longer the reads are that I produce (like simulating improvement in the long-read sequencing technology), far less SVs gets detected, while my hypothesis is that longer reads should be able to detect more. Then I found that the alignment files look clean (sam/bam), but with far more supplementary tags. I thought I could test this by changing the flags.
I am searching the internet, but changing read flags is not something they seem to do, only things like omitting reads with certain flags from the file. Maybe I am looking with the wrong search terms, in that case I am sorry.
I guess I need to change column 1 and 2 of the sam file (QNAME and FLAG). But I was also wondering if this is going to work at all. Maybe someone with more knowledge of alignment files could help me with this problem? suggestions how to change the sam/bam file to check if supplementary reads are ignored?
I think this, because I made on purpose a very unrealistic alignment. A 15Mb region against a 15Mb region with differences (x10 times to " simulate" 10 reads / 10x coverage). The alignment file looks quite good actually in IGV, but clicking on any read reveals that there are of course MANY supplementary reads. I also have a "normal" alignment file, where I also want to check if changing the supplementary tag could possibly help me with this tool.


text:
Block wise protein imputation
I am currently working on a dataset that contains 50 samples (10 samples * 5 blocks). The features of the date set are:

The data is perfectly balanced between blocks, with equal treatment representation in each block. Each block contains 2 control samples (CTL) that are pools of all other treatment samples. These are technical replicates that are identical across multiplexes.

There is a lot of missing data (~40%),

most of the missing data is missing by block, meaning that the protein is present in some blocks, but not others.

When looking at the expression of these missing proteins in the other blocks, they tend to be around the median value for the experiment, which suggests to me that these values are missing completely at random.

I could proceed using complete cases, but then I'm missing out on a large portion of data. I am considering using an imputation method, but I'm concerned that (due to the nature of the missing data) that I will just be introducing a block effect into the data.
Any input on potential methods for imputation would be greatly appreciated.
Experimental design:
SampleNumber    Multiplex   Treatment
1   A   CTL
2   A   CTL
3   A   TMT1
4   A   TMT1
5   A   TMT2
6   A   TMT2
7   A   TMT3
8   A   TMT3
9   A   TMT4
10  A   TMT4
11  B   CTL
12  B   CTL
13  B   TMT1
14  B   TMT1
15  B   TMT2
16  B   TMT2
17  B   TMT3
18  B   TMT3
19  B   TMT4
20  B   TMT4
21  C   CTL
22  C   CTL
23  C   TMT1
24  C   TMT1
25  C   TMT2
26  C   TMT2
27  C   TMT3
28  C   TMT3
29  C   TMT4
30  C   TMT4
31  D   CTL
32  D   CTL
33  D   TMT1
34  D   TMT1
35  D   TMT2
36  D   TMT2
37  D   TMT3
38  D   TMT3
39  D   TMT4
40  D   TMT4
41  E   CTL
42  E   CTL
43  E   TMT1
44  E   TMT1
45  E   TMT2
46  E   TMT2
47  E   TMT3
48  E   TMT3
49  E   TMT4
50  E   TMT4

The data is in a matrix format, columns = samples and rows = individual proteins, similar to gene expression.
The controls are technical replicates, but they are pools of all other 40 samples. This is for TMT proteomics, where the pools are used by the MS to identify peptides across the multiplexes, not necessarily for multiplex correction. I've included multiplex (as well as treatment) in the linear model to identify changes associated with the treatment, while removing the effect of the multiplex. I have confirmed (with PCA and PVCA) that the multiplex effect is negligible after including it in the lm. My problem at the moment is the missing data.
I calculated the missing percentage from the data matrix by taking the number of missing values divided by the total number of potential values. he data is in a matrix format, columns = samples and rows = individual proteins, similar to gene expression. To calculate the missing values, I divided the total number of NA's by the total number of values in the matrix (number of samples * number of proteins). Missing data is common in proteomics, but normally this can be attributed to something, most often being close to the limit of detection. But that is not the case here, if a value is missing, it is missing in the whole multiplex, including the pool, where we know it should be found.
The missing values in this case are random, and this is definitely due to the MS. It will select the first ~65k peptides that it identifies to trigger an additional MS event, so if it doesn't see that first peptide, it will be missing from the entire multiplex. There are methods for imputation of data that is missing at random, but I'm concerned that it will just be imputing by multiplex.

40% missing data is huge. Missing data analysis is complicated on the underlying distribution. If the data set is periodic then missing data periodicity is needed. Non-periodic data can be solved using regression, I think SciPy has an automated method for this, but the volume of missing data makes this approach complex. If you can bring the missing data down to 20% by compartmenatlising the data set that would be the starting point.

This looks like a TMT proteomics experiment. What tools/pipeline did you use to arrive at your quantitative values?
There are a variety of tools that are useful for analyzing this type of data and supply appropriate imputation methods. If you are comfortable programming in R, one I would suggest is the MSstatsTMT R package: http://msstats.org/msstatstmt/
That being said, any ways you can improve the upstream analysis to decrease the amount of missing values will be extremely beneficial.


text:
Bedtools wrongly indicates zero coverage
I have a BAM file produced using BWA-MEM and I need to calculate the coverage of each exon in this sequencing. I tried using this command:
bedtools coverage -a exome_sorted.bed -b exome_alignment.bam -sorted
And the result was the following:
chr1    12595   12802   ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242  .       ref|DDX11L1     0       0       207     0.0000000
chr1    13163   13658   ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242  .       ref|DDX11L1     0       0       495     0.0000000
chr1    14620   15015   ref|WASH7P,ref|NR_024540,ens|ENST00000538476,ens|ENST00000438504,ens|ENST00000488147,ens|ENST00000541675,ens|ENST00000423562       .       ref|WASH7P      0       0       395     0.0000000
...

That is, zero coverage for all exons. But there are reads mapping to the exons, as I can see by opening my BAM file on IGV and looking at some exonic regions:

I have already tried sorting the BED file using sort -k1,1 -k2,2n (as recommended on the Bedtools documentation) and using bedtools sort with the -faidx option, but neither option solved my issue.
if I skip the -sorted option, the memory usage becomes way too high.
I have also checked if the chromosome names are the same in the BAM file and in the BED file, using samtools view -H. Both use chr1,chr2,...,chr22,chrX,chrY, although BAM has chrM and stuff like chrUn_gl000214 which are not on the BED.
Thanks!
EDIT:
@ATpoint pointed out that all visible reads in this region have MAPQ0. I looked at another exon that did not have this issue, but Bedtools still outputs zero coverage. IGV picture below:

If this is your output :
chr1    12595   12802   ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242  .       ref|DDX11L1     0       0       207     0.0000000

Then I guess your bedfile looks like this (the first 6 fields) :

chrom : chr1
FeatureStart : 12595
FeatureEnd : 12802
name : ref|DDX11L1,ref|NR_046018,ens|ENST00000518655,ens|ENST00000450305,ens|ENST00000456328,ens|ENST00000515242
score : .
strand : ref|DDX11L1

The issue I see here is that you are giving the wrong information to the strand field of your bedfile. It should be either +, either - or . if you don't want to consider the strand information.
Hope it helps !


text:
Manually order the set fails in UpsetR
This question was also asked on Biostars
I have tried to use UpsetR to visualise the input file, which can be found here. How is it possible to make UpsetR accept an order defined in sets=c("Nlab", "NQLD", "Ngla", "Nsyl", "Ntom", "Ntab", "Natt", "Cann", "Stub", "Slyc", "Vvin", "Atha")?
with keep.order = T

with keep.order = F

> library("UpSetR")
> orthogroups_df<- read.table("orthogroups.GeneCount.tsv",  header=T, stringsAsFactors = F)
> #All species
> selected_species <- colnames(orthogroups_df)[2:(ncol(orthogroups_df) -1)]
> selected_species
 [1] "Atha" "Cann" "NQLD" "Natt" "Ngla" "Nlab" "Nsyl" "Ntab" "Ntom" "Slyc" "Stub" "Vvin"
> head(orthogroups_df)
  Orthogroup Atha Cann NQLD Natt Ngla Nlab Nsyl Ntab Ntom Slyc Stub Vvin Total
1  OG0000000    0    0  965    0    0    3    0    0    0    0    0    0   968
2  OG0000001    0    1    3    0    0  448    0    0    0    0    0    0   452
3  OG0000002    0    1  313    0    0  120    1    0    1    0    0    0   436
4  OG0000003    0   93   15   21   46   16   33   63   36   25   39   26   413
5  OG0000004    1   42    2   34  109    6    8  154   11    9    4    0   380
6  OG0000005    0    2   61    1   34   44   91   70   43   20    1    0   367
> ncol(orthogroups_df)
[1] 14
> orthogroups_df[orthogroups_df > 0] <- 1
> # we only show intersections of interest ,
> intersections=list(list(selected_species),
+                    list("NQLD", "Ngla", "Natt", "Nlab", "Nsyl", "Ntab", "Ntom"),
+                    list("Stub", "Slyc"),
+                    list("Atha", "Vvin"),
+                    list("Ntab", "Nsyl", "Ntom"),
+                    list("Nlab", "NQLD", "Ngla"),
+                    list("Nlab", "NQLD", "Nsyl"),
+                    list("Nlab", "Ngla", "Nsyl"),
+                    list("NQLD", "Nsyl", "Ngla"))
> upset(orthogroups_df,
+       text.scale = c(1.4),
+       sets=c("Nlab", "NQLD", "Ngla", "Nsyl", "Ntom", "Ntab", "Natt", "Cann", "Stub", "Slyc", "Vvin", "Atha"),
+       order.by = "freq",
+       keep.order=T,
+       intersections = intersections,
+       sets.x.label="Total number of orthogroups",
+       mainbar.y.label = "Number of orthogroups")

It's possible that upset is using factor levels for ordering, which is the default approach for ggplot2. If this is the case, then it can be fixed by converting the variable to a factor, and using an explicit level definition:
sets <- factor(sets, levels=sets)

For your specific example:
upSets <- c("Nlab", "NQLD", "Ngla", "Nsyl", "Ntom", "Ntab", "Natt", "Cann", "Stub", "Slyc", "Vvin", "Atha")
upSets <- factor(upSets, levels=upSets)
upset(orthogroups_df,
   text.scale = c(1.4),
   sets=upSets,
   order.by = "freq",
   keep.order=T,
   intersections = intersections,
   sets.x.label="Total number of orthogroups",
   mainbar.y.label = "Number of orthogroups")


text:
Extracting sequences from FASTA beginning with common 5' end
I am trying to figure out the best way to extract sequences from a FASTA file which begin with a common 5' region of 43 nucleotides. Preferably, I would like to to allow for "fuzziness" in this region to allow for mutations or read overlaps. The idea is that any sequences that begin with this region, regardless of size, could be extracted into a new file.
Example of the FASTA file:
>1-69050-454.08
GTACGGGGAAGGACGTCAATAGTC

>2-65989-433.95
AATCTGTAGTACGACTCACTATAGC

>3-62181-408.91
GATCTGTAATACGACTCACTATAGG

>4-49959-328.53
GGGGAAGGACGTCAATAGTCACAC

And what I would like to get from the code is such:
>2-65989-433.95
AATCTGTAGTACGACTCACTATAGC

>3-62181-408.91
GATCTGTAATACGACTCACTATAGG

The FASTA file is quite large and I have tried using a couple grep and awk methods retrieve the sequences. Any help you could provide is much appreciated.

What is the common region that you are expecting to see?  Also do you want to allow for insertions and deletions or just base edits?

Here is a 'quick' python3 solution, requiring installation of BLASTN and a pip3 install biopython.  The results are output in FASTA format to stdout.
Save the following python3 code in hit_finder.py
#!/usr/bin/env python3

import sys
from Bio import SeqIO
from Bio.Blast.Applications import NcbiblastnCommandline

cline = NcbiblastnCommandline(query=sys.argv[1],
                              subject=sys.argv[2],
                              outfmt="'6 qseqid length qstart'",
                              max_hsps=1,
                              word_size=4)

result = [j[0] for j in [i.split('\t') for i in
                         cline()[0].rstrip().split('\n')]
          if int(j[1]) >= int(sys.argv[3]) and int(j[2]) <= int(sys.argv[4])]

for record in SeqIO.parse(sys.argv[1], 'fasta'):
    if record.id in result:
        print(record.format('fasta'))

Make it executable: chmod +x hit_finder.py.
This script takes four arguments (which in future versions could be enhanced using argparse instead of sys.argv):
arg1 = query_fasta_file (file to search in)
arg2 = subject_fasta_file (database file to contain the sequence to match at 5')
arg3 = min_alignment_length (this means the length of the hit, allows mismatching in string with blastn and a word_size of 4)
arg4 = start_of_hit_at_5prime (where should the alignment hit start?)

Let's say you wanted to search for AATCTGTAGTACGACTCACTATAGC in your query seqeuences.  Save this sequence in subject.fa as:
>subject
AATCTGTAGTACGACTCACTATAGC

Your query sequences are stored in queries.fa as:
>1-69050-454.08
GTACGGGGAAGGACGTCAATAGTC

>2-65989-433.95
AATCTGTAGTACGACTCACTATAGC

>3-62181-408.91
GATCTGTAATACGACTCACTATAGG

>4-49959-328.53
GGGGAAGGACGTCAATAGTCACAC

Now, to search for hits of at least 23 bp starting no further than 1 base from the 5' end of your sequence, you would run this search using:
./hit_finder.py queries.fa subject.fa 23 1

To search for hits of at least 30 bp no more than 10 bp from the 5' end of your queries, do:
./hit_finder.py queries.fa subject.fa 30 10

To save your results, output stdout to file, e.g.,:
./hit_finder.py queries.fa subject.fa 30 10 > results.fa

This is exactly what I was looking for, I appreciate your assistance and thank you for taking your time to write this out.


text:
Generating basic dna sequences in R
Looking for some guidance.
In an introductory bioinformatics class I am tasked with doing the following in R:

Create the sequence ACTGACTG
Generate a random string of nucleotides with equal probabilities.

Are there certain functions that I should be looking at? No library or package or set of functions was provided so I am completely lost. I am assuming there must be functions in R that do these things that I need to manipulate, but I cannot seem to find any resources explaining that. I looked at seq() (thinking it might generate DNA sequences), but it only seems to generate sequences of integers. I also looked at randDNA(), but it is random, how could I make it generate the exact string in question 1).
Any direction will be greatly appreciated.
Thanks.

They are two separate problems.  If you have something that generates sequences of random integers can you think of a way to map those to the bases you want to report.  The first question is ambiguous.  How do your report the sequence you created?

Without giving away too much for your homework check out the sample() function in R. Think about how you can pass DNA bases to this function to generate random sequences.

Thanks this is really helpful. Can you please help me understand the prob parameter, as I understand the probability of each element is relative to the entire list. Is it possible to specify probabilities for each element in certain positions? say I have sample(1:2, size=4, replace=TRUE, prob=c(0.8, 0.2)) This says that 1 has probability 0.8 is the list overall, but I want to say 1 has probability 0.8 in position 4 (and some other probability in other positions)

This can't be done using just the prob parameter, because it only takes a single vector as input. You'll have to find some loopy, or otherwize crazy way for trying to apply different probabilities at different positions.


text:
Calling GEOquerry error while loading the library
Note: this question has also been asked on Biostars
I started getting a strange error when trying to load GEOquerry with library(GEOquerry).
Error: package or namespace load failed for ‘GEOquery’ in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]): namespace ‘rlang’ 0.4.9 is already loaded, but >= 0.4.10 is required

When I run > devtools::install_github('r-lib/rlang') I get ...
Error in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]) :
  namespace ‘rlang’ 0.4.9 is already loaded, but >= 0.4.10 is required

Has anyone dealt with this before?  Restarting my R session has no effect.


text:
Calculating read average length in a Fastq file with bioawk/awk
I found here this awk script:
BEGIN {
  headertype="";
  }
{
  if($0 ~ "^@") {
    countread++;
    headertype="@";
  }
  else if($0 ~ "^+") {
    headertype="+";
  }
  else if(headertype="@") { # This is a nuc sequence
    len=length($0);
    if (len>4) {
      readlength[len]++;
    }
  }
}
END {
  for (i in readlength){
    countstored+=readlength[i];
    lensum+=readlength[i]*i;
    print i, readlength[i];
  }
  print "reads read = "countread > "/dev/stderr";
  print "reads stored = "countstored > "/dev/stderr";
  print "read average length = "lensum/countstored > "/dev/stderr";
}

and I just wonder if it is possible to shorten it with bioawk?

It might be possible.  What problem do you have with having a long script? Why do you want a shorter script? Or do you want a faster script?

Did anyone get the correct script for finding out the average read length?

This can also be done with regular awk.
awk '{if(NR%4==2) {count++; bases += length} } END{print bases/count}' <fastq_file>

The NR%4==2 count the second line out of every block of 4.  length is a built-in that defaults to the length of the line, same as length($0).  In this case you can inject you custom printing to the END{} block but countread and countstored will always be the same since the averaging is done on the fly.

Note that, like OP’s code, this will only work on FASTQ files that don’t break the sequence and/or quality string across multiple lines.

That is correct.  However that practice is highly frowned upon so should not be an issue with most FASTQ files

I've yet to ever encounter a FASTQ file which has sequence or quality on multiple lines, I think for all intensive purposes assuming 4 lines per record should be fine.

This script is wrong because a quality string may start with @. With bioawk, it can be simplified to:
bioawk -c fastx '{ readlength[length($seq)]++; countread++ } END{...}'

The END{} block is the same as your original version.
EDIT: forgot -c fastx is the original answer. Thank @MattBashton for pointing this out.

Using the exact same END{} code from the OP fails in this example with bioawk: illegal field $(), name "seq"


text:
How will Seurat handle pre-normalized and pre-scaled data?
I want to use data from the following dataset:
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE84371
The data has been TPM normalized, which is not ideal for clustering but I have to work with what I have. TPM is a bad normalization method and it should not be used for these analyses because its laden with a lot of assumptions.
Presumably it has already been scaled, can I simply skip these steps and then run FindVariableFeatures() followed by the rest of the pipeline from the clustering tutorial or do the methods used assume Seurat specific preprocessing and using TPM will screw them up?
Unfortunately in this case TPM is not reversable.  https://btep.ccr.cancer.gov/question/faq/what-is-the-difference-between-rpkm-fpkm-and-tpm/#:~:text=Here's%20how%20you%20calculate%20TPM,divide%20this%20number%20by%201%2C000%2C000.
TPM is calculated as a function of the kilobases of the gene in the reference genome they used and the reads per gene across all samples (aka all cells).  They have not specified which reference genome was used in the main article, and they do not supply the read counts in the metadata or the main file.  So in this case the TPM is irreversible.
I think I can get raw reads from SRA and do alignment myself. It's annoying that raw reads were not provided but as someone who's lab entirely does computational analysis on open datasets I can tell you that this is not the most egregious thing I've ever seen. As I said before it would be easier for me to get the fastq files and do my own alignment which is still a ton of work the main authors could have prevented by simply posting the read counts.
Clustering Tutorial:
https://satijalab.org/seurat/articles/pbmc3k_tutorial.html

TPM means transcripts per million. It is obtained by dividing the read counts by the length of each gene in kilobases, resulting in reads per kilobase (RPK). All the RPK values in a sample are counted and then divided by 1,000,000.

I don't do transcriptome analysis, it ain't my thing, however I do understand statistical analysis as well as the underlying issue regarding the public availability of molecular data ... I agree with the OP its not ideal. However, yes the OP can continue with 'clustering', personally I definitely prefer it if the transformation was removed, or analysed from the raw data.
Summary The authors of the data set are obliged to provide the raw data, presuming its published. Why? This is because whilst the OP can perform the analysis - it is very likely the OP will be asked to repeat the analysis on the untransformed data (again note this ain't my area, so 'very likely' might be 'possibly' in this field). At this point the original authors have to declare their untransformed data, if this was a formal request to the OP's analysis.
Please read the caveats and discussion.

Firstly could the OP reverse the normalisation? Its just a transformation applied to all the data which in their analysis resulted in a normal distribution, I don't know TPM but if this was log-transform you would supply 'un-log' it, square-root then square it, reciprocal do another reciprocal will 'de-transformation'. Can TPM be undone? I agree it ain't ideal, from several perspectives but it gets things moving quickly.

Ask the authors for the raw data - this is what I would do. They are obliged to provide it if the data is published. Virtually all journals have this policy in place, however I agree that ASM remain an outlier in this respect (who are significant publisher). Any other journal will return to the authors and point out their 'check-boxes' when they submitted the manuscript. It is etiquette to go to the authors first and provide sufficient chance for them to respond, prior step 3 (below) - giving them the absolute maximum benefit of the doubt.

What you do if the authors don't respond? The OP would then contact the journal explaining: a) why you need the data, i.e. the analysis performed; b) explain the analytical significance of the raw data, i.e. arguments the OP has raise 'not ideal' and general concerns raised here. I do have experience in this regard and the authors will oblige if asked by the journal.

Searaut's cluster analysis certainly used tSNE ... I believe it uses dimensionality reduction such principle component analysis (PCA) - I don't know which they use - and then tSNE analysis (which I think you are calling 'clustering*'). This is a very robust 'clustering'* (please see footnote), there is no question this is as good as it gets. If that is singly your goal I agree that dimensionality reduction calculations don't need transformations, transformations are needed for parametric statistics and other areas of model fitting.

The problems with transformation
I agree the authors should not have done this because:

Transformations are not perfect, if its wrong - you're stuck.
Not only do you not need to transform the data, its better you don't - if untransformed data results in nice 'clusters' thats the preferred result. Dimensionality reduction + tSNE gives great results thats well known.

Corroborating the OP's concern
Dimensionality reduction calculations do not need a normal distribution because its sectioning up the observable variance to maximise the classification of the data. My opinion - in statical analysis - is if you need to transform, THEN perform dimensionality reduction (PCA-stuff) THEN tSNE, that data has been heavily manipulated. It doesn't mean to say the conclusions are not valid but its been put through a lot and would STRONGLY SUGGEST THE RAW DATA HAS NOT GIVEN GOOD RESULTS. Thus, personally I would be cautious, why should be OP be put under unnecessary caution?
It is important to note this is not my field and the field will have its own standards. If transformation is routine for all statical calculations (except non-parametric where it may not be valid) then authors of the data will argue they are correct and no further data needs to be released. Personally, I really doubt it, but we'll see what others say.
Footnote * clustering is a statistical method usually resulting in a tree, it is generally avoided because there are so many different clustering methods giving different results and no agreement on which is the 'optimal' method. I think when the OP say 'clustering' they refer to clusters following a tSNE analysis, which will be a bivariate plot.

From the comments Failure to declare a reference genome is a formal omission that has to be rectified and the journal is obliged to do this. Thus the OP would simply contact the journal, they may tell the OP to contact the authors and when the authors respond the response is passed back to the journal. If its an online journal you can simply leave a comment in the 'comments' section: the authors will answer. In this instance officiating via the journal is preferred to ensure the information is verified by all parties - I presume the OP does not know the authors. From the comments I get the gist the OP is working via supervision, and would therefore just notify  their supervisor who should officiate on OP's behalf (they might then delegate it back to the OP - but at least the supervisor is formally aware of this).
In context, stuff happens and there is formal process to ensure parity is achieved its simply patience in working through it.

Down vote In response to the downvote ... unpopularity happens in data set transparency its part of the process ... I was (at least in part) attributed to changing journal group policy on data availability, stated from within the journal. Maybe they joking - but they certainly all knew about what was going on. The specific issue was a prolonged dialogue over access to a specific data set. It does demonstrate that most journals are concerned if an author does not cooperate with data transparency, in fact the journal in question explicitly said this and they were not joking.

With regards to your comments about the dimensionality redcution this is par for the course in RNA-seq one point to note though is that the results of the clustering should be robust across different embeddings (tSNE, UMAP, spectral clustering, CCA, PCA etc).

Thanks for the update. Regarding the last point ... the authors of the data don't have to be cooperative, but they have to be transparent. What you have said is there is a chance they have not been transparent and that is rarely accepted with molecular data.

Seurat doesn't use tSNE for clustering; tSNE is used only for visualisation.

I notice that the paper points to the Broad website for data accession; have you looked there?
https://singlecell.broadinstitute.org/single_cell/study/SCP1/-single-nucleus-rna-seq-of-cell-diversity-in-the-adult-mouse-hippocampus-snuc-seq
For specifics on the normalisation method, you could look at their associated github repository:
https://github.com/yinqingl/nucseq_analysis
Given that the paper was published in 2016, they are TPM counts, and you have access to the raw reads, my recommendation would be to regenerate gene count tables from those raw reads. I expect it will be more time consuming, less accurate, and require more effort to track down the gene count table than to create one from scratch. You can use Salmon Alevin to do this fairly quickly:
https://salmon.readthedocs.io/en/latest/alevin.html
[bearing in mind that there may be some footwork required to work out how to get Salmon to interpret the barcoding method that was used]


text:
How to filter Ensembl cDNA and ncRNA FASTA files by primary assembly?
I'm currently performing differential expression analysis using alignment-free quantification using Kallisto. To do this, I need to create a Kallisto index using Ensembl's cDNA and ncRNA annotations available at the following two links:
ftp://ftp.ensembl.org/pub/current_fasta/mus_musculus/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz ftp://ftp.ensembl.org/pub/current_fasta/mus_musculus/ncrna/Mus_musculus.GRCm38.ncrna.fa.gz
However, later on in my analysis I realized that I had multiple Ensembl gene IDs pointing to the same gene. For example, I was gettign results for both ENSMUSG00000026012 and ENSMUSG00000102412 which both map to Cd28. I realized that the reason for this error was due to the fact that the annotations I was using contained haplotypes.
I would like to remove these haplotypes, and to do this I believe that I need to rerun my analysis using a primary assembly. However, Ensembl doesn't have annotated cDNA and ncRNA FASTA files of the primary assembly. How can I filter these files to contain only genes present in the primary assembly? Is there a better way for me to solve this issue?

Both Ensembl and GenBank databases have their own ideas about what constitutes a gen, and the mapping of Ensembl genes to RefSeq genes is not one-to-one. There are some RefSeq genes that map to multiple Ensembl genes, some Ensembl genes that map to multiple RefSeq genes, and some genes that are not in the other database.

The simplest method is to download the GTF file for GRCm38 and filter that. You can then use one of the many tools out there (bedtools getfasta (make sure to enforce strandedness and possibly use the --split option, if you happen to use a BED12 file), gtf_to_fasta from tophat2 (it's the only reason to keep tophat2 installed), etc.) to make a transcriptomic fasta file.
Alternatively, Ensembl BioMart can be used to directly download a transcriptomic fasta file using a filter to select only those transcripts on pre-defined chromosomes. Other filters can be applied as well. If you don't already have a GTF file downloaded and one of the previously mentioned packages installed, then this may be the simplest way to get what you want.

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


text:
Program to make a haplotype network for a specific gene
Currently my lab manually uses pop art to make the network after using PHASE to id an individual's alleles. I would like to know how to generate the graph structure (i.e. the links between vertices), as well as how to draw the graph.
Input data are formatted as FASTA files with the sequence in nucleotide formatting.
Example:

I just feel like there's got to be a better, more time effective way to do this.

This is a haplotype map, where each node is proportional to the frequency of a given allele. The formal tool to do this is Network 5 found here. The problem with manaul haplotype maps is the area of the node should be proportional to the frequency, if this is not accurate it could be perceived to be misleading.
The approach for looking at each parsimonious state along a branch is a phylogenetics method which is performed via MacClade 4.0 ... requires an old Mac though.

What do network edges represent?


text:
"Error in col2rgb(colors, alpha = TRUE) : invalid RGB specification" when running HeatmapAnnotation() (from ComplexHeatmap)
This question was also asked on github
I'm trying to run the function HeatmapAnnotation (a function in the ComplexHeatmap package). When I run it, I get the error "Error in col2rgb(colors, alpha = TRUE) : invalid RGB specification"
col2rgb() is not among the libraries I've loaded and it doesn't seem to be a method contained in HeatmapAnnotation.
Here's the code block for the section:
df <- data.frame(labid= combmetadata$Sample_Name,
             Diagnosis= combmetadata$Diagnosis,
             Tissue = combmetadata$Tissue,
             Tumor.Normal = combmetadata$T.N,
             study = combmetadata$batch)

 labid <- colnames(bValsselect)
 df <- left_join(as.data.frame(labid),df,by="labid")
 df <- df[,c(2,4,5)]
 df$study[df$study == "1" | df$study == "2"] <- "cgid"

 df$Diagnosis <- factor(df$Diagnosis, levels = c("Uterine endometrial","Ovarian endometrioid",
                                            "Ovarian mucinous",
                                            "Colorectal mucinous",
                                            "Pancreas mucinous",
                                            "Stomach mucinous"))

 ha_column = HeatmapAnnotation(df = df,
                          col = list(Diagnosis = c("Uterine endometrial" = "#DDCC77F",
                                          "Ovarian endometrioid" = "#0F7554",
                                          "Ovarian mucinous" = "#44AA99",
                                          "Colorectal mucinous" = "#882255",
                                          "Pancreas mucinous" = "#AA4499",
                                          "Stomach mucinous" = "#D695D0"),
                            Tumor.Normal = c("N" = "grey",
                                             "T" = "black"),
                            study = c("cgid" = "#332288",
                                      "tcga" = "#BD9AF7")))

hegt <- Heatmap(bValsselect,
            column_split=4,
            name = "Beta Value",
            col = colorRamp2(c(0,0.25,0.5),c("#00FFCC","#FFFFFF","#0099FF")),
            show_row_names = FALSE,
            top_annotation = ha_column,
            show_row_dend = F)

The code runs fine up until it reaches the end of the HeatmapAnnotation() call
I've tried to load col2rgb() and override the alpha setting. Can anyone give me a clue about what's happening and any suggestions on how to correct the problem?  I'm running R 4.2. Thanks

Looking at the documentation, it seems that HeatmapAnnotation expects a vector, whereas you are giving it a list. It is not clear on what acceptable color types are. It might even work with a list, but I'd suggest playing around with what you are supplying as your col parameter.
Note that HeatmapAnnotation has to do color mapping of some kind, the code is a little dense, but it does look like it calls out to SingleAnnotation, which in turn runs ColorRamp2, which does use col2rgb.
It also looks like col2rgb is called in ColorRamp2, which you use in the Heatmap call. How sure are you that your script isn't making it that far? It would help to see error output, possibly.
That said, when I take your isolated ColorRamp2 and run it, it goes ok. Possibly pull out the ColorRamp2 call and run it yourself in your environment to make sure nothing weird is happening.
Try moving the assignment of the colors outside the Heatmap call, occasionally supplying function calls as parameters like that causes issues.

You're right about the chain of calls.  I do know that the code didn't run past the HeatmapAnnotation() function because in the process of debugging it, I could generate the error before I run the Heatmap() subblock.  However, I just realized that the problem is that the first color assignment contains 7 alphanumeric characters instead of 6.  That's what was causing the function to fail.  Thanks for your help though.  It certainly helps me understand the function a little better.

The alphanumeric value for the first color assignment contains 7 alphanumeric characters instead of 6. Correcting this typo fixed the problem.


text:
Will blastn remove sequences from a search with low identity?
I'm using this command (excuse the duplicate naming, I know it's bad form):
blastn -query cm-seqs/combined_seqs.fna -db combined_seqs.fna -out cm-matched.txt -num_alignments 1 -outfmt 10

to take a set of fasta sequences to blast against another data base I built with blast.
My query file has 4,364,417 sequences. The files resulting from the above command spat out 4,362,639 results. The difference is about 1800. Is this due to poor alignments for some sequences that didn't pass any of blasts default parameters? I assumed the result would be 1 "best" match in the database for each query sequence.

-num_alignments 1 is not appropriate for that output format i think, try -max_target_seqs 1

Yes, if a query sequence has no good alignment, blast will not return a result for it. Your command will indeed output only the best hit for each query sequence, but there is no reason to assume all of your sequences will have a match.
Those that don't have a match won't be included in the output. This sounds perfectly normal unless you are blasting the same input file against a DB created from that file. Otherwise, you don't expect all of your sequences to necessarily have a match.
Also note that since the default blastn settings are quite permissive (usually an e-value of 10, for example), don't make the mistake of assuming that a blast hit means anything. The details will always depend on the specific search you are running (I have discussed this a little in my answer here), but in the vast majority of cases, any hit with an e-value > 1 can safely be ignored. In most cases, you are only interested in e-values well under 1, actually. So you probably already have all sorts of spurious results in there anyway.
For more details, please edit your question and explain what the query and DB files are.

Great! That was a perfect response.

The databases I created were the bacterial 16S metagenome of an additive, and the 16S metagenomic community of our experimental subject during the experiment. The microbiome of the subjects were blasted against the additive to track what in the additive was able to persist in a foreign environment.

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


text:
Samtools/bcftools calling indels in noisy reads
I have very noisy nanopore reads and am trying to call SNPs/indels.
I'm runinig into some trouble when truing to use the samtools mpileup | bcftools call combination.
It seems to incorrectly be calling very long indels where it looks like there is no support. For example:
samtools mpileup -f ref.fa sample.bam -r Chromosome:198940-198940 produces:
Chromosome    198940    C    37    .,.-1A.-1A...-1A.-1A.-1A.,-1a.-1A.,-1a,-1a.,.-1A,,.-1A,,+1a,+1a,.-1A,,-1a.,-1a,-1a,,-1a,-1a.-1A,-1a.    06?5/<:8>/4<?3691465<3354.08:23:11@3C

This looks fine... however when I try to call the genotype with bcftools I get a very long indel.
samtools mpileup -ugf ref.fa sample.bam -r Chromosome:198940-198940 -t AD | bcftools call -mv produces:
Chromosome    198940    .    CAAAAGACGTCATCGACGTACGGTCGGTGACCACCGAGATCAACACGTTGTTCCAGACGCTCACCTCGATCGCCGA    C    225    .    INDEL;IDV=32;IMF=0.5;DP=62;VDB=0.145385;SGB=-0.662043;MQSB=1;MQ0F=0;AC=2;AN=2;DP4=0,0,7,2;MQ=60    GT:PL:AD    1/1:255,27,0:0,9

Using htsbox pileup -f ref.fa sample.bam | grep 198940 it shows there is one read with this indels together with many other reads with smaller indels.
Chromosome  198940  C   C,C-75AAAAGACGTCATCGACGTACGGTCGGTGACCACCGAGATCAACACGTTGTTCCAGACGCTCACCTCGATCGCCGA,C-4AAAA,C-2AA,C-1A,T,C+1A,C+1T,C+3GGT 0/4:25,1,1,2,24,5,2,1,1

Is there any reason why it would be grouping these reads together into the big indel?
Any idea what is going wrong here?

It's likely that the caller is being overly sensitive to deletions, which are the most common error mode for nanopore reads, even three years later in 2021.
What I'd guess is happening is that the caller is seeing a string of consecutive 1bp deletions at high frequency, and treating these as a consensus deletion covering many more bases. Looking at the pileup calls for the surrounding bases may help to diagnose this.
My general advice to people doing variant calling with nanopore is to be suspicious of deletions, especially if they include homopolymers. That usually means trying to confirm any important variants with another technology.
Newer, more accurate, base calling models are being created all the time by ONT, so if you still have the raw FAST5 data, it might be worth it to run through the most recent base caller.


text:
How are the values of prop.part() and prop.clades() calculated?
Consider the following dataset:
fictional.df <- data.frame(L1 = c(0,0,0,0,0,0,0,0),
                       L2 = c(0,1,0,0,0,1,1,0),
                       L3 = c(1,1,0,1,1,1,1,1),
                       L4=c(0,0,1,1,0,0,0,0))

I

converted this to a phyDat object and then
created a pairwise distance matrix as follows:

    fictional.phydat <- as.phyDat(fictional.df,
                              type="USER",levels=c("1","0"),
                              names=names(fictional.df))
    fictional.hamming <- dist.hamming(fictional.phydat)

From this distance matrix, I then estimated a UPGMA tree:

    fictional.upgma <- upgma(fictional.hamming)

I then created bootstrap datasets:

    set.seed(187)
    fictional.upgma.bs <- bootstrap.phyDat(fictional.phydat, FUN =
    function(xx) upgma(dist.hamming(xx)), bs=100)

I then calculated the proportion of partitions in the bootstrap set:

    upgma.bs.part <- prop.part(fictional.upgma.bs)

So far so good. Here is where I would appreciate some help. When I call the function prop.clades, I do not understand the result:

    prop.clades(fictional.upgma,fictional.upgma.bs)
    [1] 100  NA  71

Question Why does this function return NA when there is evidence for that clade in the set of bootstrap trees?
A second question:
prop.clades(fictional.upgma,part=upgma.bs.part)
[1] 100  49 112

If there are only 100 bootstrap samples, why is the value for the final clade 112?


text:
Is there an R package that computes homoplasy excess ratios (HER)?
I want to calculate the homoplasy excess ratio for a phylogenetic tree, for example homoplasy resulting from recombination.
I am wondering if there are any R packages out there that will do this?
I have looked at phyloseq, Rphylopars, phylocurve, PhylogeneticEM, PhyloinformR, but these packages do not have functions that calculate HER. I searched "homoplasy excess ratio" on bioconductor, but it returned ca. 3000 packages.

Maynard's Smith's homoplasy is what you want and it is available in the recombination detection package (RDP) described here. Its really a Windows package rather than R.