# Question
Title: Is there a way to get a list of all PDB IDs for proteins with crystal structure information?
I would like to download all the PDB files which have crystal structure information in them. But the query builder seems to require me to enter the PDB IDs manually.

Is there a way to get a list of those PDB IDs at RCSB PDB?

Thanks in advance.

# Answer
> Comment migrated as answer as it inadvertently may have been an answer.

If the sole records in a PDB file that are of interest are `EXPDTA`,`CRYST1` ORIGX1\`\`,`ORIGX2`, `ORIGX3`,`SCALE1`,`SCALE2`,`SCALE3`, the in the FTP site the file ftp.wwpdb.org/pub/pdb/derived\_data/index/crystal.idx contains this information in a tabular format.

If an extra piece of metadata is sought, such as **some of** the refinement details in `REMARK` lines, then the PDBe API can help. In particular, the `experiment` entrypoint (example) is of interest.
> 3 votes

# Answer
running the query builder at RCS PDB against id with space group (so crystal info) I get:

more about RCSB PDB web services here: https://www.rcsb.org/docs/programmatic-access/web-services-overview

and the search API here: https://search.rcsb.org/index.html#search-api think you should query for entries where `symmetry.pdbx_full_space_group_name_H_M` exists (https://search.rcsb.org/structure-search-attributes.html and ask to have returned entry.id (a list of in this case)
> 2 votes

---
Tags: pdb
---


# Question
Title: Search ChebiEntity search by common or synonym name?
I can download information from ChEBI if I know the corresponding ChEBI ID:

```
from libchebipy import ChebiEntity
chebi_entity = ChebiEntity("15903")

```

And then look at all the synonyms that the compound has:

```
chebi_entity.get_names()
> [{'_Name__name': 'beta-D-Glucose', '_Name__typ': 'SYNONYM', '_Name__source': 'KEGG COMPOUND', '_Name__adapted': False, '_Name__language': 'en'},
 {'_Name__name': 'beta-D-glucose', '_Name__typ': 'NAME', '_Name__source': 'UniProt', '_Name__adapted': False, '_Name__language': 'en'},
 {'_Name__name': 'beta-D-glucopyranose', '_Name__typ': 'IUPAC NAME', '_Name__source': 'IUPAC', '_Name__adapted': False, '_Name__language': 'en'},
 {'_Name__name': 'BETA-D-GLUCOSE', '_Name__typ': 'SYNONYM', '_Name__source': 'PDBeChem', '_Name__adapted': False, '_Name__language': 'en'},
 {'_Name__name': 'WURCS=2.0/1,1,0/[a2122h-1b_1-5]/1/', '_Name__typ': 'SYNONYM', '_Name__source': 'GlyTouCan', '_Name__adapted': False, '_Name__language': 'en'}]

```

Is there a way to get the ChEBI ID with one of the synonyms?

```
beta-D-glucose -> 15903

```

Is there any other python tool or web-api that can achieve this?

# Answer
> 1 votes
How very interesting ...

Yes its easy just download the flat file database here:

> https://ftp.ebi.ac.uk/pub/databases/chebi/Flat\_file\_tab\_delimited/names.tsv.gz

**WARNING**: it's a 5.5M file ... don't click the link unless you really want to download a file this size.

You can simply load it straight into a `pandas` dataframe from `read_csv()` and then interrogation it via `loc`. You could place the header `NAME` as the index. There are other way in, you can shift the dataframe to a dictionary and interrogate that way, or place it straight it SQL via e.g. `sqlite3`

Here I just loaded it onto Excel and found the entry you were looking for ...

....

---
Tags: embl
---


# Question
Title: Retrieve rsIDs from chromosome positions using MRutils::get_rsid_from_position
I have 5609 chromosome positions with both ref and alt alleles. I want to retrieve rsIDs for the same. I tried `get_rsid_from_position()` of the `MRutils` package in R, but after some time while fetching data from the internet, the cache memory of my system got full, and it froze the system. Any other way to resolve this issue?

I have 200 GB allotted for the cache memory.

# Answer
> 3 votes
An alternative possibility in R is using the library biomaRt from Bioconductor. Here an example starting from a table with 3 columns: Chr, Start, End. Organism is human and positions are assumed to be based on the reference version GRCh38.

```
library(biomaRt)
snpMart = useEnsembl(biomart = "snps",
         dataset = "hsapiens_snp")

# assuming your table is tab-separated, with column names
myData <- read.table("myTable.txt", sep="\t", header=True)

coords <- apply(myData, 1, paste, collapse = ":")

# submit the query
getBM(attributes = c('refsnp_id', 'chr_name', 'chrom_start', 'chrom_end', 'allele'),
filters = c('chromosomal_region'),
values = coords,
mart = snpMart)

```

Source

---
Tags: r, snp, rsid
---


# Question
Title: How to subset an SRA file for a single chromosome?
I used prefetch to get the Pacbio reads of chicken from the SRA database. I want to align these reads against a reference genome, but not all the reads. I am only interested in a particular region on Chromosome Z. Is there a way that I can subset the downloaded SRA file?

I am trying to use blast, but don't know where to start. How should I prepare the SRA file from prefetch for the blast? I have currently started fasterq-dump for the downloaded SRA file.

Or should I subset after aligning the entire SRA reads data with the reference genome, and in that case what tools do I use?

# Answer
> 3 votes
After dumping out the reads using `fasterq-dump`, you'll need to align them first and then extract those that map to your region of interest. I think minimap2 is an excellent choice. If you have HiFi reads, you can just use something like:

```
minimap2 -ax map-hifi ref.fa pacbio-ccs.fq.gz > aln.sam

```

Then, using samtools, coordinate-sort and index the sorted alignments:

```
samtools sort -o aln.sorted.bam aln.sam
samtools index aln.sorted.bam

```

With an indexed BAM, you can then view the alignments that overlap your region of interest:

```
samtools view aln.sorted.bam chrZ:1-10000

```

If you instead write these to a file, you can then convert them using samtools fasta. I think something like the following should suffice:

```
samtools view -bo chrZ_1_10000.bam aln.sorted.bam chrZ:1-10000
samtools fasta -0 chrZ_1_10000.fa -o /dev/null chrZ_1_10000.bam

```

You could then upload the FASTA file to the NCBI Web BLAST server.

---
Tags: assembly, reference-genome, pacbio, sratoolkit
---


# Question
Title: Examples of machine learning approaches to validate results where ground truth is lacking?
I am currently looking into some of the published methods for deconvolution of spatial transcriptomics data where each spot does not have single-cell resolution. These methods all rely on cell type annotated scRNA-seq data to build a reference of expected expression for each cell population. So in essence, because ground truth is more or less impossible to obtain for complete cell typing of a spatially intact tissue section, another type of data (scRNA) is used instead. This lets you infer something about the spatial organization of the tissue despite no spatial information being present in the training data, but in the end no ground truth is available for validating the results.

I was trying to think of an analogous machine learning problem, biological or not, where something similar is done to deal with unattainable ground truch. Are there any obvious/good examples of something like this? Maybe some inspiration can be drawn for how such results are validated.

# Answer
> 1 votes
If the problem concerns spatial clustering of cells then the ML algorithm is without question DBSCAN.

The problem with deep learning and CNN (convolution) not "deconvolution", i.e. the MBE paper cited, is:

* The model is an *ad hoc* construction, which can be also a strength;
* Under-training, is resolved by adding more layers;
* Overtraining, can only be assessed post-analysis, in ML its algorithm specific;
* CNN approaches require large datasets for training of 2 million+.

The issue with ML is its performance tapers at around the 1 or 2 million data set size mark, whereas deep learning *can* (stress *can*) continue to use massive data sets to enhance performance (thats its strength).

Whether, DBSCAN has received publications in bioinformatics - I'm not sure. However, what is present in the public domain isn't always the cutting edge work: a lesson I've very recently learnt. DBSCAN is very fashionable for spatial problems at present and in context personally would rank it above CNN.

**CNN vs ML** The thing is most of us - unless you're working for Google - don't have access to data sets in 2 million+ range. This can be circumvented by *transfer learning* but it is not sufficiently developed outside Google and a few major players in computer vision. Has Google performed CNN on one million+ cell images - I doubt it. If the problem is ecological species identification via remote cameras - Googles CNN image libraries would be just the ticket (lots of dogs, cats and bear photos) but for microscopy studies I seriously doubt it.

**Conclusion** What I am 100% certain on is DBSCAN is VERY fashionable right now and certain ML algorithms repeatedly outperform deep learning in certain situations. What I am trying to say is ML is not necessarily just a low data size version of deep learning: in some instances it surpasses it and I think DBSCAN has that potential. The broader challenges of deep learning hinge around 'transfer learning'. If biology manages to resolve that *then* deep learning will start surpassing machine learning, but firstly it's about sharing, secondly in e.g. RNN this is not technically trivial.

---

What you can do is pose a detailed problem in a later question and I could explain how training would be used ... generally in complex analytical problems it is difficult to assess precisely how a solution would be implemented.

---
Tags: rna-seq, scrnaseq
---


# Question
Title: Improving prokaryotic assembly with other contig/scaffold-level data?
I have what at first sight appears to be a high-quality MAG (~10 pieces, high completion%) that I built from a hybrid assembly (Illumina + Nanopore data) from a cyanobacterium.

Workflow: Quality control (BBDuk) \> Hybrid assembly (Unicycler) \> co-binning multiple different samples (vamb) \> Quality control (CheckM, anvi'o)

\*The Nanopore data is pretty old, and therefore I'd thought Unicycler might be a safer bet rather than long-read-first assemblers, due to the higher error rate

The same strain also has 1-2 very fragmented genomes (100s of pieces) available in NCBI at scaffold level. The reads themselves are not available, so I can't just add them to the assembly.

Is there a way to use the available NCBI genome(s) to try to check/polish mine? (I've checked out Improve a reference genome with sequencing data but it's already a few years old so I imagine different tools might be applicable)

Thank you for your time!

# Answer
> 2 votes
If your goal is to compare to the NCBI version I'd suggest using QUAST to compare the two assemblies, supplying the NCBI assembly as the reference genome (`-r` option I think). This will not guarantee anything, but you would expect the two assemblies to be fairly syntenic and to have good overlap of sequences. The tool includes visualization, but you could also use an approach like circoletto or dot plots for something a little more intuitive/simpler.

I'd expect your genome to be quite a bit higher quality even as a metagenome-derived assembly due to nanopore contiguity, I'm not sure that using the NCBI reference to improve yours is a great idea. I'd instead suggest following @M\__'s suggestion to scaffold with your nanopore reads (you will have to map them onto your full metagenome first, though approach may differ based on what the scaffolder wants). I'm not up on the current methods here but here is a recent tool. If you have 10 pieces I'd expect to be able to close at least some of that using nanopore (even old).

If you do choose to use the reference to improve, the tool that people seem to be using now is RagTag.

# Answer
> 2 votes
**Summary** What I am saying in this post is there approximates to 3 "generations" of *in silico* gap-closers/polishes:

1. GapFiller - amongst other "early" generations (2012);
2. medaka/racon/Nanopolish/Pilon (2014 - 2018 inclusively)
3. Nextpolish/PEPPER/Apollo/Homopolish/NeuralPolish (2020 onwards).

The trendy approach is combining 2. and 3. I can't guarantee it will work for cyanobacteria gap-closing, but the shear volume of algorithms dedicated to this purpose, infers they must be delivering results.

**Background** Whats the issue with a reference genome? The key thing about *many* bacteria is their ability to "dump" genes and then regain them via plasmid associated recombination. This can occur extensively within a single species. This makes reference associated assembly essentially impossible in many bacteria - *particularly against the official NCBI reference genome*. Hence, the official *NCBI reference gene*, is not like a human genome reference - which is a gold standard. With due to respect to those who have dedicated time to presenting the *official bacterial reference genomes* are not very useful.

This BioStars post here demonstrates that even *Mycobacteria tuberculosis* would not assemble via a reference genome and TB in context is one of the most conserved genomes in bacteria, the SNP behaviour is famously limited.

If attempting a reference genome it is worth considering the following:

1. Any amount of genetic divergence can result *genetic island* phenomena and non-contiguous genomes that will make assembly difficult, therefore a complete genome as close to the target genome as possible is a good reference.

* For example, simply Blasting sections of you genome to look for a consistent reference that is close to yours. A custom database of complete genomes for your species my be useful, but just sending bits of your target onto the NCBI website will be useful. The idea is to get lots of bits of the genome matching the closest strain.

2. QUAST was going to be my suggestion too.
3. Mapping via e.g. Mummer4 against the reference genome as the "scaffold".

**Key question:**

> The same strain also has 1-2 very fragmented genomes (100s of pieces) available in NCBI at scaffold level. The reads themselves are not available, so I can't just add them to the assembly.

I'm not sure what "not available" means.

**Key suggestion: Polishing** The other approach are recent *de novo* assemblers and *in silico* "polishing". What I'm recommending is to first try "polishing". This doesn't exclude subsequent reference genome assemblies - *not at all* \- but it's a way to improve the *de novo* assemblies and if there are still issues then consider a reference genome.

Firstly I would read someone who is *really* up to date, e.g.

> Comparative evaluation of Nanopore polishing tools for microbial genome assembly and polishing strategies for downstream analysis

They are medaka fans and are using this as part of a combination approach. medaka was part of the established stable of "polishers"

* GapFiller \- I don't think this is useful.
* racon is specifically for long reads
* Nanopolish here

The authors prefer medaka then goes on to describe

> Nextpolish, PEPPER, Apollo, Homopolish, and NeuralPolish

Much more recent from 2020 ... thats a lot of polishing tools. Note, PEPPER and Apollo are often used for algorithms names - make sure its 2020 onwards, with PEPPER being 2021. They suggest medaka x PEPPER.

I really think this combination polishing approach with an established (medaka) and one the very latest polishing algorithms, such as PEPPER (2021) is really worth trying. If it doesn't work for the specific purpose it will certainly improve the contig quality. However, given the shear volume of work dedicated to this area then there has to be good reason.

By the same token that

> Money doesnt give happiness, but sure gives a better standard of misery.

Then ...

> Polishing might not give gap closure, but sure gives a better standard of quality to the same old bunch of miserable contigs

---

***de novo* Assembly** Maybe Miniasm? I'd look at "polishing" the *de novo* assembly you have.

RagTag as a SPAdes alternative is interesting @MaximilianPress

---
Tags: genome, assembly, nanopore, illumina
---


# Question
Title: How to get separate histograms plots on the basis of the column value?How to detect which plot has most deviation?
I have a data frame (df) which has correlations calculated for different genes with respect to different ID combinations. I want to get separate histogram plot based on the gene name (separate plot for separate gene). Also after plotting I want to know which genes have more deviations in their correlation values (Like more deviated from the median in each plot).

Example: `df`

```
Gene    ID_1    ID_2    Correlation
TNF 1175077258  1175075956  0.626115
TNF 1175077261  1175075956  0.542499
TNF 1175077261  1175077258  0.50113
TNF 1175077263  1175075956  0.431587
TNF 1175077263  1175077258  0.734858
TNF 1175077263  1175077261  0.622206
CCL2    952399885   952399878   0.632703
CCL2    952399894   952399878   0.622244
CCL2    952399894   952399885   0.697313
CCL2    952399923   952399878   0.687826
CCL2    952399923   952399885   0.612089
CCL2    952399923   952399894   0.603834
IL10    1068768000  1068767927  0.118955
IL10    1068768147  1068767927  0.32038
IL10    1068768147  1068768000  0.430287
IL10    1068768409  1068767927  0.335264
IL10    1068768409  1068768000  0.406426
IL10    1068768409  1068768147  0.546452

```

I tried the below code which only works when combining all genes together.

```
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams["figure.figsize"] = (8,6)

g1 = sns.histplot(df.Correlation)
g1.set(
    xlim=(-1,1),
    title="correlation"
)
g1.axvline(0.0, c='black', linestyle="--")
g1.axvline(df.Correlation.median(), c='red', linestyle="--")

```

The result/plot I am expecting is this (I have created this example plot in excel)

# Answer
> 2 votes
You could use the groupby() pandas function to group the dataframe by gene name. And then just loop through each group to plot the histogram of correlation values. I think you can use the mean of the ratio (deviation / median) to compare deviation for different genes:

```
import matplotlib.pyplot as plt

for name, group in df.groupby('Gene'):

    median = group['Correlation'].median()
    deviation = abs(group['Correlation'] - median)
    ratio = deviation / median
    ratio_mean = ratio.mean()

    id_1 = group['ID_1'].astype(str)
    id_2 = group['ID_2'].astype(str)

    group['Combination'] = id_1 + '/' + id_2

    plt.figure(figsize=(8, 6))
    plt.bar(group['Combination'], group['Correlation'])

    plt.xlabel('Combination', fontsize=12, fontweight='bold')
    plt.ylabel('Correlation', fontsize=12, fontweight='bold')

    plt.xticks(rotation=45, ha="right")
    plt.title(f'{name} (deviation ratio: {ratio_mean.round(3)})')
    plt.tight_layout()

    plt.gca().spines['top'].set_visible(False)
    plt.gca().spines['right'].set_visible(False)

    plt.show()

```

---
Tags: python, visualization, correlation, pandas
---


# Question
Title: How to get a file with the introns of a gene?
I was able to get the exons of a particular gene as a fasta file from the NCBI genome viewer. As for introns I only found a bed file that showed where the intron starts and ends, along with its direction. What tool do I use to get the exact sequence of these 11 introns using this bed file? And is there a more direct way? Thank you.

# Answer
First I edited the Chromosome column of the bed file from ncbi, because it contained the accession number instead of the chromosome. Then I used the getfasta tool of the bedtools utility.

`bedtools getfasta -fi input_genome_file.fa -fo output.fa -bed bed_file`

The output contained the introns in the format I wanted. If anybody had a better way, feel free to post.
> 1 votes

---
Tags: gene, reference-genome, bed
---


# Question
Title: How to find specific types of assemblies for specific species using entrez tools?
How to find specific types of assemblies for specific species using entrez tools?

Task: Trying to specifically find transcriptomes and associated cDNA data for a list of speices.

I can use this oneliner `esearch -db assembly -query "{species name }" | esummary` to get a list of assemblies and associated data available for any given species.

I can use `esearch -db assembly -query "{species name }" | xtract -pattern ENTREZ_DIRECT -element Count` to find any specific info nested within the XML tree.

Now I am not sure how to find if the assembly is a transcriptome (which is what I am looking for). I have tried to `CRTL +f` terms like `rna`, `peptide`,`amino acids` to no avail and I guess there is a chance I am missing something while looking it over manually.

For example here is the output for "zea mays".

Does anyone have any clues? Or should I be using a different db or different tool?

# Answer
I strongly recommend switching to SRA https://www.ncbi.nlm.nih.gov/sra

Whether the following would work ...

```
esearch -db sra -query "{species name }" | esummary

```

I'm not at all sure because SRA has its own toolkit. Maybe NCBI `e-utilities` `e-summary` `efetch` would work on the SRA database. However, the weblink certainly works for transcriptomics.

NCBI has shifted some of its database to be housed off-site onto GCP and AWS and I have a feeling these may not be accessible to `e-utilities`, but I could be wrong.

Just to note, what I envisage the question meaning is a complete transcriptomic study for a given species (or given sample) under given experimental conditions.

---

Comments. I better clarify, SRA will be raw Illumina output, it is not assembled into exome cDNA. `e-utils` is basically for assembled cDNA, it might give you the information - I really don't think it will download it though. So to answer your questions, yes the data will be there on SRA, whether this is the format you want is a different question.
> 0 votes

---
Tags: assembly, transcriptome, entrez
---


# Question
Title: How to statistically test the strength and weakness of a correlated pair?
Subset of `df` before `groupby` looks like

```
Diabetes_Group  DM-1    DM-2    corr
Group1  56  56  1
Group1  56  58  0.611
Group1  56  61  0.448
Group1  56  63  0.3587
Group1  58  56  0.611
Group1  58  58  1
Group1  58  61  0.60113
Group1  58  63  0.534858
Group1  61  56  0.448
Group1  61  58  0.60113
Group1  61  61  1
Group1  61  63  0.622206
Group1  63  56  0.3587
Group1  63  58  0.534858
Group1  63  61  0.622206
Group1  63  63  1
Group2  78  78  1
Group2  78  85  0.622703
Group2  78  94  0.622244
Group2  78  23  0.687826
Group2  85  78  0.622703
Group2  85  85  1
Group2  85  94  0.697313
Group2  85  23  0.612089
Group2  94  78  0.622244
Group2  94  85  0.697313
Group2  94  94  1
Group2  94  23  0.603834
Group2  23  78  0.687826
Group2  23  85  0.612089
Group2  23  94  0.603834
Group2  23  23  1
Group3  27  27  1
Group3  27  80  0.118955
Group3  27  147 0.32038
Group3  27  9   0.335264
Group3  80  27  0.118955
Group3  80  80  1
Group3  80  147 0.430287
Group3  80  9   0.406426
Group3  147 27  0.32038
Group3  147 80  0.430287
Group3  147 147 1
Group3  147 9   0.546452
Group3  9   27  0.335264
Group3  9   80  0.406426
Group3  9   147 0.546452
Group3  9   9   1

```
1. What statistical test can be used to do comparative analysis between the different ID's (`DM1 or DM2`) to check which ID pair has weaker/stronger correlation within each diabetes\_group?
2. Is there any statistical test to categorise the data into three categories:

(a) `pairs of DM1 AND DM2` which are highly correlated than others.

(b) pairs that are not highly correlated

(c) For the category 1 pairs (highly correlated), is there any particular ID that outperforms than other ID's among different groups?

# Answer
This comprehensive analytical outline for the stats and all the stats questions could be captured inside 3 analyses. Implementing this will obviously take a bit of time, but there should not be any hidden caveats (transformations etc ...).

**Question 3**

The statistical test you want is K-means clustering based on the correlation score where the number of predefined groups is 3. The labels ... you will need to combine the group with the

```
Label          Value
Group1-56-56   1
...

```

This is an absolute solution to the data and really cool use of K-means.

**Question 3c** An easy solution is:

* perform a new K-means with e.g. 10 groups and find the group with highest correlation values.

A better solution is:

* hierarchical clustering would identify specific groups and for this sub-question is a good analysis. The input is the correlation matrix, the output is a tree. HOWEVER, interpreting this could be difficult if you're not trained in reading trees.

---

**Question 2**

A 2-way ANOVA would be a really good solution to this data for Q2. The test would be the columns were DM1 DM2 and the rows Group1, Group2, Group3. Interpreting it would be fun because you'd also have to grapple with a concept called "interaction" - which could be important here. However, the ability to address a whole raft of questions within *one* analysis is what 2-way ANOVA would permit. There is a slight issue of unequal group sizes 1, 2 and 3, because it would be better if Group 1,2 and 3 etc... were equal in size. However, the raw data looks to be very amenable to ANOVA. You must not use the correlation index in this test - just the raw values.

---

K-means clustering is in scipy and `scikit-learn`

Hierarchical clustering is also in scipy

For ANOVA

```
import statsmodels.api as sm
from statsmodels.formula.api import ols

```

There will be easy approaches to 2-way ANOVA e.g. SPSS and I think Excel (needs a plug in). I suspect that Hierarchical clustering will also be in a nice point and click package like SPSS (I don't use them, so I don't precisely know.)
> 2 votes

---
Tags: python, statistics
---


# Question
Title: Can I traverse/enumerate PDB one at a time using a Python script without downloading all of them to my local disk?
What would be the total size of all PDB files in RCSB?

Can I traverse/enumerate them one at a time using a Python script without downloading all of them to my local disk?

If YES, what would be the Python script like?

# Answer
> 6 votes
The page https://www.wwpdb.org/ftp/pdb-ftp-sites gives the addresses of interest. PDB has an FTP site downloadable via rsync and the latter has a dry run mode.

```
rsync --port=33444 -rlpt -v -z --dry-run --stats rsync.rcsb.org::ftp_data/structures/divided/pdb/

```

You get a bunch of folders names and the following:

```
Number of files: 182479
Number of files transferred: 0
Total file size: 38382243431 bytes
Total transferred file size: 0 bytes
Literal data: 0 bytes
Matched data: 0 bytes
File list size: 3777395
File list generation time: 52.909 seconds
File list transfer time: 0.000 seconds
Total bytes sent: 16
Total bytes received: 3777651

sent 16 bytes  received 3777651 bytes  66861.36 bytes/sec
total size is 38382243431  speedup is 10160.30

```

So there are 182,479 \*.ent.gz in various folders. The total is 38,382,243,431 bytes, which is 38 GB, so not too much. In Python it is easy to iterate across gzip files with the gzip module, so uncompressing them is generally a thing to avoid as it will bloat up greatly. The problem arise when dealing with the density maps, which do take up a lot of space...

## Download without saving

The PDBe is easier to for downloading things than the RCSB PDB and the API are much more comprehensive.

Given a list of PDB codes one could

```
pdbcode : str = ''

import requests
r = requests.get(f'https://www.ebi.ac.uk/pdbe/entry-files/download/pdb{pdbcode.lower()}.ent', stream=True)
assert r.status_code == 200, f'something ({r.status_code} is wrong: {r.text}'
pdb_code = r.text

```

# Answer
> 0 votes
Not Python, but you can do this with BioStructures.jl in Julia:

```
using BioStructures
pdblist = pdbentrylist()
for pdbid in pdblist
    # Download a file and delete it after use
    downloadpdb(pdbid, format=MMCIF) do fp
        s = read(fp, MMCIF)
        # Do something with the file, e.g. print the cysteine fraction
        nres = countresidues(s, standardselector)
        if nres > 0
            frac = countresidues(s, standardselector, x -> resname(x) == "CYS") / nres
            println(pdbid, "  ", round(frac; digits=2))
        end
    end
end

```

If you want to get all the PDB entries you should use the mmCIF format, though it does lead to larger file sizes.

# Answer
> 0 votes
my attempt working as of the 17th January 2023 using RCSB PDB Search API:

```
import json

# import urllib.parse

import requests

import sys

# data = {"query":{"type":"terminal","service":"text","parameters":{"attribute":"rcsb_entry_container_identifiers.entry_id",
            # "operator":"exists","negation":False}},"return_type":"entry"}

data = {"query":{"type":"terminal","service":"text","parameters":{"attribute":"rcsb_entry_container_identifiers.entry_id",
             "operator":"exists","negation":False}},"request_options":{"return_all_hits":True},"return_type":"entry"}

data = json.dumps(data)

headers = {'content-type': 'application/json'}
url = 'https://search.rcsb.org/rcsbsearch/v2/query?json='

params = {}

risposta = requests.post(url, params=params, data = data, headers=headers)

# print(risposta.status_code,'\n', risposta.content.decode())

rispo = json.loads(risposta.content.decode())

# print(rispo["result_set"], type(rispo["result_set"]), len(rispo["result_set"]))

id_list = []
for i in rispo["result_set"] :

    id_list.append(i['identifier'])

print(id_list, type(id_list) , 'lenght : ', len(id_list), 'size bytes  : ' ,sys.getsizeof(id_list))

```

you get a huge list :

```
.............................
'2V0H', '2V0G', '2V0F', '2V0E', '2V0D', '370D', '2ULL', '2UKD', '2UGI', '2UDP', '2UCZ', '2UBP', '2UAG', '2TYS', '369D', '368D', '367D', '366D', '2TUN', '2TSS', '2TSR', '2TSY', '365D', '2TRT', '2TRS', '2TS1', '2TRY', '2TRX', '364D', '2TSC', '2TSB', '2TSA', '363D', '2TRC', '2TRA', '2TRM', '2TRH']

<class 'list'> lenght :  200069 size bytes  :  1671784

```

Using the requests GET endpoints :

```
import requests

import sys

url = 'https://search.rcsb.org/rcsbsearch/v2/query?json=%7B%22query%22%3A+%7B%22type%22%3A+%22terminal%22%2C+%22service%22%3A+%22text%22%2C+%22parameters%22%3A+%7B%22attribute%22%3A+%22rcsb_entry_container_identifiers.entry_id%22%2C+%22operator%22%3A+%22exists%22%2C+%22negation%22%3A+false%7D%7D%2C+%22request_options%22%3A+%7B%22return_all_hits%22%3A+true%7D%2C+%22return_type%22%3A+%22entry%22%7D'

risposta = requests.get(url)

print('+++++++++++++++++++++++++++++')

rispo = risposta.json()

id_list = []
for i in rispo["result_set"] :

    id_list.append(i['identifier'])

print(id_list, type(id_list) , 'lenght : ', len(id_list), 'size bytes  : ' ,sys.getsizeof(id_list))

print('--------------------------------------')

print(risposta.status_code)

print('\n\n',risposta.url)

```

getting as of the 19th January 2023 --\> 200375 entries

---
Tags: python, pdb, biopython
---


# Question
Title: Velvet Optimizer automatically changes to hash-length 31
I'm trying to use Velvet Optimizer for a De Novo Assembly; I set my hash-lengths to be between 55 and 69

```
/projects/path/to/perlscripts/VelvetOptimiser-2.2.6/VelvetOptimiser.pl -s 55 -e 69 -f '-shortPaired -fastq.gz -separate trial2.trimmed.paired.R1.fastq.gz trial2.trimmed.paired.R2.fastq.gz' -f '-short -fastq.gz trial2.trimmed.unpaired.R1.fastq.gz' -f '-short -fastq.gz trial2.trimmed.unpaired.R2.fastq.gz' -x 2

```

When I check the output file, it returns me with

```
Starting to check input parameters.
    Start hash value too high.  New start hash value is 31.
    End hash value not in workable range.  New end hash value is 31

```

Any idea where I could be going wrong in my command?

# Answer
If a value for `MAXKMERLENGTH` was not set when compiling `velveth`, VelvetOptimizer just uses the default `$maxhash` value of 31. The hash lengths that you specify on the command line must be less than or equal to this value to run the tool successfully. The solution is to re-compile `velveth` with a larger `MAXKMERLENGTH`, for example:

```
make 'MAXKMERLENGTH=69'

```

Alternatively, re-install using the package manager:

```
conda install -c bioconda perl-velvetoptimiser

```

This package includes velvet, which is built using a much higher maximum k-mer length and is usually sufficient for most purposes: `MAXKMERLENGTH=191`.
> 4 votes

---
Tags: assembly
---


# Question
Title: Is it possible to do homology inference across species using different kinds of NGS data?
Background: I have a list of species that I want to put through homology inference. The goal of homology inference is to investigate the evolution of a trait on a species tree. I want to use the quantitative aspect of the data like CNV and recreate nodal CNV using cafe or count. And use gene trees of the homology to see if I can discern divergence patterns across gene families of interest.

The issue is some species only have RNA-seq generated transcriptions (peptides and cDNA) while others have genomes. The issue of missing introns aside, how can I circumvent this gap?

Should I:

1. Just use the cDNA data with genomic data? compare nucleic acid to nucleic acid. I need to eventually make gene trees as well and as such MSA is important.The list of species is huge and I anticipate seg-faults during MSA if not during homology inference itself.
2. Use emboss or something else to convert the nucleic acid data to peptides.And then compare peptides to peptides? Which makes the homology inference → MSA (hopefully pal2nal and gene tree inference later) a little easier. I have never done translation before so I do not know the pitfalls per se.

My list of species is green plants.

Can anyone share their experience if you have done something like this before?

# Answer
> 4 votes
Personally, I would use `blastx` against a local protein database. This converts query cDNA-\>amino acids and therefore combines points 1 and 2 without need to perform a translation in advance.

Genetic distance and saturation might become an issue here depending whether the search is outside the family and especially Order

---

MAFFT vs Muscle Super5. I would try "Muscle5" or rather the latest Muscle under the "Super5" algorithm. It outperforms everything.

---
Tags: rna-seq, phylogenetics, phylogeny, sequence-homology
---


title: dealing with a list of VEP files

text:
# Question
Title: dealing with a list of VEP files
I'm new to dealing with big and compressed data so, I have some questions regarding that

I have 1000 files, each file has the extension (vep.txt.gz) I want to read these files, then filter them. The size of these files is very large. I tried to read these files in Python and R but I am facing a memory error issue. I need to read these files and extract some rows from each file for example I have two columns. In the first column, I want the rows with low values in the impact column and the rows with likely pathogenic values from the clinSIG column. Also, I have to add a new column to put the files' names as IDs for each sample.

So, how can I read a list of big (vep.txt.gz) files in R or Python, then filter the content of the files based on two columns (impact = low) and (clin\_SIG = Likely pathogenic) ?

Should I decompress the files or just read them?

The following is an example of the files but each file is compressed and all the files have the same header (you can find the original file, which is tab separated, \[here\]\[1\]):

```
## ENSEMBL VARIANT EFFECT PREDICTOR v108.0
## Output produced at 2017-03-21 14:51:27
## Connected to homo_sapiens_core_108_38 on ensembldb.ensembl.org
## Using cache in /homes/user/.vep/homo_sapiens/108_GRCh38
## Using API version 108, DB version 108
## polyphen version 2.2.2
## sift version sift5.2.2
## COSMIC version 78
## ESP version 20141103
## gencode version GENCODE 25
## genebuild version 2014-07
## HGMD-PUBLIC version 20162
## regbuild version 16
## assembly version GRCh38.p7
## ClinVar version 201610
## dbSNP version 147
## Column descriptions:
## Uploaded_variation : Identifier of uploaded variant
## Location : Location of variant in standard coordinate format (chr:start or chr:start-end)
## Allele : The variant allele used to calculate the consequence
## Gene : Stable ID of affected gene
## Feature : Stable ID of feature
## Feature_type : Type of feature - Transcript, RegulatoryFeature or MotifFeature
## Consequence : Consequence type
## cDNA_position : Relative position of base pair in cDNA sequence
## CDS_position : Relative position of base pair in coding sequence
## Protein_position : Relative position of amino acid in protein
## Amino_acids : Reference and variant amino acids
## Codons : Reference and variant codon sequence
## Existing_variation : Identifier(s) of co-located known variants
## Extra column keys:
## IMPACT : Subjective impact classification of consequence type
## DISTANCE : Shortest distance from variant to transcript
## STRAND : Strand of the feature (1/-1)
## FLAGS : Transcript quality flags

#Uploaded_variation  Location   Allele  Gene             Feature          Feature_type  Consequence                            cDNA_position  CDS_position  Protein_position  Amino_acids  Codons   Existing_variation  IMPACT    DISTANCE  STRAND  CLINSIG
11_224088_C/A        11:224088  A       ENSG00000142082  ENST00000525319  Transcript    missense_variant                       742            716           239               S/I          aGc/aTc  -                   LOW  -         -1      likely pathogenic
11_224088_C/A        11:224088  A       ENSG00000142082  ENST00000534381  Transcript    downstream_gene_variant                -              -             -                 -            -        -                   MODIFIER  1674      -1      -
11_224088_C/A        11:224088  A       ENSG00000142082  ENST00000529055  Transcript    downstream_gene_variant                -              -             -                 -            -        -                   HIGH  134       -1      -
11_224585_G/A        11:224585  A       ENSG00000142082  ENST00000529937  Transcript    intron_variant,NMD_transcript_variant  -              -             -                 -            -        -                   MODIFIER  -         -1      -

```

From the previous file, the expected output will be one row which is the first row because the impact is low & CLINSIG is likely pathogenic so i want to find all rows with this condition for 1000 files

```
ID  #Uploaded_variation  Location   Allele  Gene    Feature Feature_type  Consequence   cDNA_position   CDS_position    Protein_position  Amino_acids  Codons   Existing_variation  IMPACT    DISTANCE  STRAND  CLINSIG
M1457.vep   11_224088_C/A   11:224088   A   ENSG00000142082 ENST00000525319 Transcript  missense_variant    742 716 239 S/I aGc/aTc -   LOW  -  -1  likely pathogenic

```

Thanks in advance

# Answer
> 3 votes
For these requirements, it might be possible to simply use filter\_vep included in the VEP package for this. Note that multiple `--filter` flags are ANDed together (i.e. all filters must pass for a line to be printed). For example:

```
for i in *.vep.txt.gz; do
    filter_vep \
        -i "${i}" \
        -o "${i%.vep.txt.gz}.filtered.vep.txt" \
        --filter "IMPACT is LOW" \
        --filter "CLIN_SIG matches likely_pathogenic"
done

```

If you then wanted to merge the output and include the filename, you could just use . With the following in a file called `merge_filtered_vep_txt.awk`:

```
BEGIN {
    FS=OFS="\t"
    f=1
}

/^##/ {
    next
}

sub(/^#/, "") && f {
    print "ID", $0
    f=0
    next
}

{
    print FILENAME, $0
}

```

Run using:

```
awk -f merge_filtered_vep_txt.awk /path/to/*.filtered.vep.txt

```

Python can also open gzip-compressed files for reading (and writing), using the gzip module in the standard library. These can be parsed line-by-line to avoid memory errors. Here's one way using the csv module (also in the standard library):

```
import argparse
import contextlib
import csv
import gzip
import pathlib
import sys

```

```
def get_argument_parser():

    parser = argparse.ArgumentParser(add_help=False)

    group = parser.add_argument_group('input options')
    group.add_argument(
        "-i",
        "--input",
        dest='glob_pattern',
        type=str,
        default="**/*.vep.txt.gz",
        metavar="STR",
        help="The relative glob pattern (default: '%(default)s')",
    )
    group.add_argument(
        "-d",
        "--dir",
        dest='glob_path',
        type=pathlib.Path,
        metavar="DIR",
        default='.',
        help="Glob files in DIR (default: '%(default)s')",
    )

    group = parser.add_argument_group('output options')
    group.add_argument(
        "-o",
        "--out",
        type=pathlib.Path,
        metavar="FILE",
        default='-',
        help="Write the filtered output to FILE (default: stdout)",
    )
    group.add_argument(
        "-f",
        "--force",
        action='store_true',
        help="Overwrite the output file if it exists",
    )

    group = parser.add_argument_group('additional options')
    group.add_argument(
        "-h",
        "--help",
        action="help",
        help="Show this help message and exit",
    )

    return parser

```

```
def main():

    parser = get_argument_parser()
    args = parser.parse_args()

    with contextlib.ExitStack() as stack:
        if args.out.name != '-':
            out = stack.enter_context(open(args.out, 'wt' if args.force else 'xt'))
        else:
            out = sys.stdout

        # Get a list of files from the directory specified by args.glob_path

        file_list = sorted(args.glob_path.glob(args.glob_pattern))

        # Build up a unique list of header fields from each input file

        header_cols = ['ID']
        for gz_file in file_list:
            with gzip.open(gz_file, 'rt') as infile:
                for line in infile:
                    if not line.startswith('##'):
                        break

                assert line.startswith('#')

                for field in line[1:].rstrip().split('\t'):
                    if field not in header_cols:
                        header_cols.append(field)

        # Create a csv.DictWriter() and write the header fields

        writer = csv.DictWriter(out, fieldnames=header_cols, delimiter='\t')
        writer.writeheader()

        # Now re-open each file in the arguments list for reading

        for gz_file in file_list:
            with gzip.open(gz_file, 'rt') as infile:

                # Skip the info header lines
                line = next(infile)
                while line.startswith('##'):
                    line = next(infile)

                # Get the column headers
                fieldnames = line.rstrip().removeprefix('#').split('\t')

                # Parse each input file and store the desired rows in memory
                reader = csv.DictReader(infile, fieldnames=fieldnames, delimiter='\t')
                rows = []
                for row in reader:

                    impact = row.get('IMPACT')
                    clinsig = row.get('CLINSIG')

                    if impact == 'LOW' and clinsig == 'likely pathogenic':
                        rows.append({'ID': gz_file, **row})

                # Write the list of rows to the output file
                writer.writerows(rows)

```

```
if __name__ == '__main__':
    main()

```

With the above in a file called `filter_vep_txt.py`:

```
usage: filter_vep_txt.py [-i STR] [-d DIR] [-o FILE] [-f] [-h]

input options:
  -i STR, --input STR  The relative glob pattern (default: '**/*.vep.txt.gz')
  -d DIR, --dir DIR    Glob files in DIR (default: '.')

output options:
  -o FILE, --out FILE  Write the filtered output to FILE (default: stdout)
  -f, --force          Overwrite the output file if it exists

additional options:
  -h, --help           Show this help message and exit

```

# Answer
> 1 votes
You don't need R or Python for this, it is a very straightforward text parsing problem and you can easily do it in an awk one-liner:

```
$ awk 'BEGIN{FS=OFS="\t"}
      /^[^#]/ && $14=="LOW" && $17="likely pathogenic" ' file
11_224088_C/A   11:224088   A   ENSG00000142082 ENST00000525319 Transcript  missense_variant    742 716 239 S/I aGc/aTc -   LOW -   -1  likely pathogenic

```

**Explanation**

* `BEGIN{FS=OFS="\t"}`: the `BEGIN{}` block is run before processing the file. Here we are setting the field separator (`FS`) and output field separator (`OFS`) to a tab so that awk will split each input line on tab characters.
* `/^[^#]/ && $14=="LOW" && $17="likely pathogenic"`: this will evaluate to true if the current line does not start with a `#`, and the 14th field is `LOW` and the 17th field is `likely pathogenic`. In awk, the default action when something evaluates to true is to print the current line, so this results in the expected output.

Now, if you also want to add the file name as a field, you can do this to have it as the last field:

```
$ awk 'BEGIN{FS=OFS="\t"} /^[^#]/ && $14=="LOW" && $17="likely pathogenic"{$(NF+1)=FILENAME; print} ' file
11_224088_C/A   11:224088   A   ENSG00000142082 ENST00000525319 Transcript  missense_variant    742 716 239 S/I aGc/aTc -   LOW -   -1  likely pathogenic       yourfilename

```

Next, if your files are compressed and you want to do this to all of them, you could simply do:

```
for file in *vep.txt.gz; do
  zcat "$file" | awk -v fname="${file%.vep.txt.gz}" '
                       BEGIN{FS=OFS="\t"}
                       /^[^#]/ && $14=="LOW" && $17="likely pathogenic"{
                          $(NF+1)=FILENAME;
                          print
                       }'
done > allfiles.filtered

```

---
Tags: r, python, filtering, vep
---


# Question
Title: How to z-transform Fst and -log(p) values for genome wide selection scan?
I saw in the literature a lot genome wide scans showing $F\_{ST}$ values in a 'z-transformed' format which seems to be bounded between 0 to infinity.

Examples are found easily searching for 'genome-wide scan $F\_{ST}$' and here is a paper where I could find such z-transformed $F\_{ST}$:

1. Identification of genomic variants putatively targeted by selection during dog domestication (see figure 1)
2. Genomic divergence during feralization reveals both conserved and distinct mechanisms of parallel weediness evolution (see figure 3)

How are the z-transformed $F\_{ST}$ values calculated?

The only reference I could find so far is this: Calculate $F\_{ST}$to draw Manhattan plot. And this Population analysis of the Korean native duck using whole-genome sequencing data where they just take the ($F\_{ST}$ \- mean $F\_{ST}$)/sd\_$F\_{ST}$. But this is giving a lot of negative values, so should this be also surrounded by an absolute value?

Then is it possible to use the z-scores to get -log(p-values)?

# Answer
> 3 votes
The standard score can be estimated by using the sample mean and sample standard deviation as estimates of the population values. Just as you wrote it in the question.

The z-score is given by

$ z = \frac{x - \bar{x}}{s}$

where:

$x$ is the sample observation
$\bar{x}$ is the mean of the sample,
$s$ is the standard deviation of the sample.

Negative values mean that particular observation is $z$ standard deviations below the mean value, and negative $z$-scores are very much expected. The process of standardization is a simple linear transformation and thus does not remove or change any degree of skewness in the original data. I guess your data is skew from the begging to values below the mean.

Yes you can compute the exact $p$-values from the $z$-scores and then apply a negative log transform. You can use a $z$ pre-computed table for this. Check this video for an example.

# Answer
> 1 votes
@JRodrigoF is the correct answer. I'll attempt to explain why the z-score is being used, because it is not immediately obvious.

$F\_{ST}$ is not trivial for anyone to understand *per se*. It is extremely useful for assessing migration - but what is the actual index that Ron Fisher devised?

Simulations suggest that $F\_{ST}$ = 0.2 is threshold above which which strong population structure is considered. Thus $F\_{ST}$ = 0 panmixia and $F\_{ST}$ = 1 is zero migration. How to interpret everything in-between is difficult, except for the 0.2 threshold if the theoretical assumptions in the threshold calculation apply to real-world data. Thus the Z-score is an interesting concept.

In other words, the 0.2 threshold simulation suggests $F\_{ST}$ is a non-linear index, which is very interesting in understanding it. Non-linear indices are a bit hairy because we all think they are linear, so 0.1 to 0.3 is one increment of migration, 0.3 to 0.5 another increment of migration. This is not true according to the simulation studies.

Z-scores get around this difficult by simply standardising the combined $F\_{ST}$ and thats why they are using it.

It becomes particularly interesting if $F\_{ST}$ is calculated per allele because then "fitness traits" can then be tracked. The allele doesn't necessarily need to be the actual fitness trait because of molecular "hitchhiking" (old John Maynard Smith terminology). Thus $F\_{ST}$ Z-score within a population within a genome is the beginning of a datamining process. So example, the domestication of dogs why it becomes such a big deal, particularly if the comparison is a dog breed that is closest to the wolf genome. The actually traits selected for in dogs can now be identified. It would also operate between genomes, to identify outliers in dog dispersal, but the behavioural (behavioral) modifications that make dogs dogs is very interesting.

---

The question is are super small probabilities generated to make log meaningful? I'm not sure in this instance probabilities well beyond very highly, highly significant are generated. If they are then a log transformation is useful.

For background the log-scale is used for transforming very small probabilities rather than use *e* nomenclature. It makes them more understandable for comparative purposes and allows you to do cool things like add and subtract them instead of multiplying and dividing. The smaller the probability value the higher the transformed value ... so in this case the higher the value the more interesting the result.

However, if the field simply does not use log transformation then it's not cool, because no-one will be able to comparatively interpret the transformed values. In phylogenetics it's done all the time, because the probabilities are really small, but anyone can immediately understand the values - because everyone uses this system. A key part of stats is clarity, not necessarily for those outside the field, but certainly for those within the field.

---
Tags: statistics, gwas
---


# Question
Title: How to fix the Error in Kruskal Wallis rank sum test for comparing the guides in python? Also to check whether using the right test for my data?
`df`

```
gene_name   guide_1 guide_1_new Correlation
MMP-1   A   A   1
MMP-1   A   B   0.426115
MMP-1   A   C   0.522499
MMP-1   A   D   0.431587
MMP-1   B   A   0.426115
MMP-1   B   B   1
MMP-1   B   C   0.60113
MMP-1   B   D   0.534858
MMP-1   C   A   0.522499
MMP-1   C   B   0.60113
MMP-1   C   C   1
MMP-1   C   D   0.622206
MMP-1   D   A   0.431587
MMP-1   D   B   0.534858
MMP-1   D   C   0.622206
MMP-1   D   D   1

```

I used four guides in an experiment (data = `df`). I calculated correlation between them based on the their log fold changes. I want to check which guide is outperforming or underperforming in comparison to other three guides in my experiment. I tried Kruskal and Wallis test to do this analysis with the script given below. However, I got an error `ValueError: Samples must be one-dimensional`. Also is this the right test to do this analysis?

`My script`

```
from scipy import stats

x = df[['guide_1']]
y = df[['guide_1_new']]
z = df[['Correlation']]
stats.kruskal(x, y, z)

```

# Answer
Reformatting the data for the KW and doing the test is easy.

Using "guide\_1" alone and pivot() to get A = \[1, 0.426115 ...

```
df2 = df.pivot(columns='guide_1', values='Correlation')

```

This will now present A, B, C, D etc ... into columns...

```
x = df2[['A']].dropna()
y = df2[['B']].dropna()
z = df2[['C']].dropna()
a = df2[['D']].dropna()

stats.kruskal(x, y, z, a)

```

That will work ... whether that is a good test is a complicated question because this is log transformed data. Transformation is used in parametric tests. Removing the log transformation is strongly recommended because KW works with raw data.

You can switch the calculation to `guide_1_new` instead. It will give a different answer ... but I don't know the experiment. The original guide looks cool.

Looking at the data - the test should work because I believe KW is 'unpaired' (you need to check) and looking at the value the distribution looks pretty similar between A, B, C and D. However, just because the KW says 'yes' doesn't mean thats a robust answer (log data). If the untransformed data says 'yes' ... thats a robust answer for KW. Just check the 'unpaired' stuff - thats important for the result.

---

Personal note, I originally used a third df but was not needed,

```
df3 = pd.concat([df2[col].sort_values().reset_index(drop=True) for col in df2], axis=1, ignore_index=True)

```
> 3 votes

---
Tags: statistics, pandas
---


# Question
Title: To get background_gradient plot for all data frames in the for loop
I tried the below code to get the background\_gradient plot for each diabetes\_group. It is giving me only plot for the last group in the plot. How to get background\_gradiant for all items (`x`) in `diabetes_group` with item name on the top ?

Subset of `df` before `groupby` looks like

```
Diabetes_Group  DM-1    DM-2    corr
Group1  56  56  1
Group1  56  58  0.611
Group1  56  61  0.448
Group1  56  63  0.3587
Group1  58  56  0.611
Group1  58  58  1
Group1  58  61  0.60113
Group1  58  63  0.534858
Group1  61  56  0.448
Group1  61  58  0.60113
Group1  61  61  1
Group1  61  63  0.622206
Group1  63  56  0.3587
Group1  63  58  0.534858
Group1  63  61  0.622206
Group1  63  63  1
Group2  78  78  1
Group2  78  85  0.622703
Group2  78  94  0.622244
Group2  78  23  0.687826
Group2  85  78  0.622703
Group2  85  85  1
Group2  85  94  0.697313
Group2  85  23  0.612089
Group2  94  78  0.622244
Group2  94  85  0.697313
Group2  94  94  1
Group2  94  23  0.603834
Group2  23  78  0.687826
Group2  23  85  0.612089
Group2  23  94  0.603834
Group2  23  23  1
Group3  27  27  1
Group3  27  80  0.118955
Group3  27  147 0.32038
Group3  27  9   0.335264
Group3  80  27  0.118955
Group3  80  80  1
Group3  80  147 0.430287
Group3  80  9   0.406426
Group3  147 27  0.32038
Group3  147 80  0.430287
Group3  147 147 1
Group3  147 9   0.546452
Group3  9   27  0.335264
Group3  9   80  0.406426
Group3  9   147 0.546452
Group3  9   9   1

```

`My code`

```
corr_matrix = []

for x, group in df.groupby('diabetes_group'):
    matrix = group[["diabetes_group", 'DM-1', 'DM-2', 'corr']]
    matrix_pivot = pd.pivot_table(matrix, values = 'corr', index=['DM-1'], columns=['DM-2'])
    d = matrix_pivot.style.background_gradient(cmap='Set3')

```

# Answer
> 1 votes
```
corr_matrix = []

for x, group in df.groupby('diabetes_group'):
    matrix = group[["diabetes_group", 'DM-1', 'DM-2', 'corr']]
    matrix_pivot = pd.pivot_table(matrix, values = 'corr', index=['DM-1'], columns=['DM-2'])
    d = matrix_pivot.style.background_gradient(cmap='Set3').set_caption(x)
    display(d)

```

# Answer
> 1 votes
Sorry I am really late answering, the `set_caption` was cool.

**Background** What we are talking about here is the `style` Class in `pandas` of which `set_background` is one method. Usually, it is for outputting a `pandas` dataframe as a webpage. If I was front-end developer it would be cool. In this example it is used to colour the different bits of the dataframe, presumably to make it easier to read a lot of data.

The thing I would have done is dump to Excel. This will preserve the colouration. The bits contained within `#` are my codes.

```
from Pathlib import Path
import pandas as pd

##########
path = Path('/path/todir/')
##########
for x, group in df.groupby('diabetes_group'):
    matrix = group[["diabetes_group", 'DM-1', 'DM-2', 'corr']]
    matrix_pivot = pd.pivot_table(matrix, values = 'corr', index=['DM-1'], columns=['DM-2'])
    d = matrix_pivot.style.background_gradient(cmap='Set3').set_caption(x)

    ##########
    pathout = Path(path, x + '.xlsx')
    d.to_excel(pathout)
    ##########

```

What was happening originally was the 'loop' was overriding the diabetes `Group` in previous loop, hence appeared to be disappearing. What will happen now is an excel file will be produced separately for each group and the file name labeled as `Group1.xlsx`, `Group2.xlsx` etc ...

Here I dumped my test output to HTML simply via `to_html` and then to pdf.

**Notes** The variable `corr_matrix = []` does not take part in this code. It may be important to other bits of code. I'm not sure `matrix` variable does much, but I've not checked.

The issue with IPython `display` is the output is not saved, so the code needs to be run each time the analyst wants to examine the correlation coefficients between `DM-1` and `DM-2`.

---

My personal dataframe was

```
df = pd.DataFrame({'diabetes_group':['Group1','Group1','Group2','Group2', 'Group3','Group3'],'DM-1':[56,56,56,56,56,56], 'DM-2':[56,58,61,63,56,58], 'corr':[0.2,0.3,0.5,0.3,0.2,0.1]})

```

---
Tags: python, visualization
---


# Question
Title: How to highlight the specific peptide sequences after performing multiple sequence alignment for the fasta file?
I have the peptide sequences and fasta files separately. I first aligned the fasta files using `msa` package. After that I'm trying to highlight the peptide sequences in the multiple sequence alignment output. I couldn't find a way to do that. Any suggestions on how to do this would be very helpful.

My R script so far:

# load required packages
library(seqinr)
library(msa)
library(dplyr)

# read the file containing the list of peptide sequences
peptides <- read.table("./peptides/AB_hydrolase-1_domain-containing_protein.txt", header=F, stringsAsFactors=F)

# read the file containing the list of fasta files
#hemoSeq <- readAAStringSet(system.file("examples/HemoglobinAA.fasta",package="msa"))
fasta_file <- "./fasta/AB_hydrolase-1_domain-containing_protein.txt"
fasta_data <- readAAStringSet(fasta_file)

# multiple sequence alignment with ClustalW algorithm
msa <- msa(fasta_data, method = "ClustalW")
```

# Answer
> 1 votes
Checkout an MSA visualisation and analysis package that you can load the alignment into. It may not be the newest tool, but I've used Jalview (https://www.jalview.org/) in the past

---
Tags: r, bioconductor, multiple-sequence-alignment
---


# Question
Title: How to rename a file based on a directory name?
I have a lot of FASTA files that I recently scaffolded with RagTag. The output of RagTag includes a directory for each sample and the corresponding FASTA file written as "ragtag.scaffold.fasta"

I would like to go through each file and rename it like so:

```
Original Path with File name:
directoryName_ID/ragtag.scaffold.fasta

Changed File name:
directoryName_ID/ID_ragtag.scaffold.fasta

```

As you can see, the directory name includes the ID number I would like to use to rename the fasta file with (along with an underscore to space it out). Is there a way to loop through all of the directories, using that last number after the "\_" in the directory name and to use that to rename my fasta file?

Here is what I have tried so far:

```
for filename in ./directory*/ragtag.scaffold.fasta; do
    mv filename ${./directory%*_*}_ragtag.scaffold.fasta;
done

```

It doesn't work at all but I know I need to use mv to rename the filename. However, when it comes to referencing the directory name the file is located in I get lost. I also have a .txt file that contains a list of IDs if that might be easier to reference. I'm still kinda new to the command line so if anyone could help that would be greatly appreciated!

# Answer
Just execute the following immediately above the directories of interest, (**post edited** to include OPs "constant term")

```
for d in *; do
  mv "$d"/ragtag.scaffold.fasta
 "$d/${d/000mergedCluster_/}_ragtag.scaffold.fasta";
done

```

The key thing I assumed is that the directory name comprises a "constant term" and is suffixed by a "variable term".

Just to note the `$d` here is directory not file name - because that is constant in each directory.

---

Comments...

```
$d/${d/000mergedCluster_/}

```

The key part of the code is above.

What its doing is saying from the variable `$d` delete the term `000mergedCluster_` (or whatever the "constant term" is)
> 3 votes

# Answer
Assuming the output is only a list of directories and each directory includes a file called "ragtag.scaffold.fasta", you could use brace expansion to generate the SOURCE and DEST strings:

```
for i in directoryName_* ; do
    echo mv "${i}"/{,"${i#directoryName_}"_}ragtag.scaffold.fasta
done

```

But the above will only work using Bash and Zsh, and maybe some other shells. Without Bash/Zsh, just use:

```
fasta="ragtag.scaffold.fasta"
for i in directoryName_* ; do
    echo mv "${i}/${fasta}" "${i}/${i#directoryName_}_${fasta}"
done

```

Remove the `echo` to perform the rename.
> 3 votes

# Answer
If you want to use Perl:

```
#!/usr/bin/perl
use 5.030;
use warnings;

use File::Find::Rule qw( );

for my $file (File::Find::Rule->name("ragtag.scaffold.fasta")->in(".")){
  my $newfile = $file;
  $newfile =~ s{(.+_(\d+)/)(.+.fasta)}{$1$2_$3};

  if (-e $newfile) {
    warn "can't rename $file to $newfile: $newfile exists\n";
  }elsif (rename $file, $newfile) {
    # do nothing
  }else{
    warn "rename $file to $newfile failed: $!\n";
  }
}

```
> 1 votes

# Answer
Perl one-liner:

```
ls */filename | perl -lane 'print "mv $1/$3 $1/$2_$3" if /(\S+_([^\/_\s]+))\/(\S+)/'

```
> 1 votes

---
Tags: fasta, bash, command-line, ragtag, for-loop
---


# Question
Title: How to write FASTA records using "Bio.SeqIO.write()"
**Update:** Biopython document says that "Bio.SeqIO.FastaIO.FastaWriter" class is obsolete. Now my question becomes how to I use `Bio.SeqIO.write()` fucntion to achive the same thing I asked.

I use the following code to extract specific FASTA records from a file based on a list of FASTA IDs given in another file. The script works just fine. However, I want to able to write the results using "Bio.SeqIO.FastaIO.FastaWriter" class into a file instead of Python's print or write function. How can I do that in my last if statement?

```
from Bio import SeqIO
#import pprint
import textwrap

my_id_file     = open('sample_files2/fasta_id_records.txt','r')
my_fasta_file  = open('sample_files2/fasta_file.fasta','r')

my_dictionary = {} # fasta IDs are keys, value can be anything.
for line in my_id_file:
    my_dictionary[line[:-1]] = 'value'
#pprint.pprint(my_dictionary)

for seq_record in SeqIO.parse(my_fasta_file, "fasta"):
    if seq_record.id in my_dictionary:
        sequence = str(seq_record.seq)
        fasta_record = textwrap.fill(sequence, width=60)
        print(f">{seq_record.id}\n{fasta_record}")

my_fasta_file.close()
my_id_file.close()

```

# Answer
> 3 votes
I figured it out:

```
from Bio import SeqIO
#import pprint
import textwrap

my_id_file     = open('sample_files2/fasta_id_records.txt','r')
my_fasta_file  = open('sample_files2/fasta_file.fasta','r')
result_file    = open('sample_files2/results.fasta', "w")

my_dictionary = {} # fasta IDs are keys, value can be anything.
for line in my_id_file:
    my_dictionary[line[:-1]] = 'value'

#pprint.pprint(my_dictionary)
for seq_record in SeqIO.parse(my_fasta_file, "fasta"):
    if seq_record.id in my_dictionary:
        #use this line to write standard fasta record using SeqIO.write
        SeqIO.write(seq_record, result_file, "fasta")

my_fasta_file.close()
my_id_file.close()
result_file.close()

```

# Answer
> 2 votes
You could instead pass in a generator that yields SeqRecords. This avoids storing the entire list of sequences in memory and is usually the way to go assuming that you need to handle large files containing thousands of records. This example use a generator expression, which is a bit like a list comprehension, but instead of creating a list, it creates a generator:

```
import argparse
import contextlib
import pathlib
import sys

from Bio import SeqIO

```

```
def get_argument_parser():

    parser = argparse.ArgumentParser()

    parser.add_argument(
        "fasta_file",
        type=pathlib.Path,
        metavar='FILE',
        help="Extract sequence records from FILE",
    )
    parser.add_argument(
        "ids_file",
        type=pathlib.Path,
        metavar='FILE',
        help="Sequence record IDs/names to extract, one per line",
    )
    parser.add_argument(
        "-o",
        "--output",
        type=pathlib.Path,
        metavar='FILE',
        default="-",
        help="Write extracted sequences to FILE (default: stdout)",
    )
    parser.add_argument(
        "-f",
        "--force",
        action='store_true',
        help="Overwrite the output file if it exists",
    )

    return parser

```

```
def main():

    parser = get_argument_parser()
    args = parser.parse_args()

    with contextlib.ExitStack() as stack:

        ids_file = stack.enter_context(open(args.ids_file))
        fasta_file = stack.enter_context(open(args.fasta_file))

        if args.output != '-':
            out = stack.enter_context(open(args.output, 'w' if args.force else 'x'))
        else:
            out = sys.stdout

        record_ids = {line.rstrip('\n') for line in ids_file}

        fasta = SeqIO.parse(fasta_file, "fasta")
        records = (rec for rec in fasta if rec.id in record_ids)

        SeqIO.write(records, out, "fasta")

if __name__ == "__main__":
    main()

```

# Answer
> 2 votes
my attempt at changing the way Biopython formats fasta file sequence width, read into code to figure out standard behaviour:

```
from Bio import SeqIO, __version__
#import pprint
# import textwrap

print('Biopython version : ' , __version__)

# my_id_file     = open('sample_files2/fasta_id_records.txt','r')
# my_fasta_file  = open('sample_files2/fasta_file.fasta','r')
#result_file    = open('sample_files2/results.fasta', "w")
### CHANGED TO:

my_id_file     = open('fasta_id_records.txt','r')
my_fasta_file  = open('fasta_file.fasta','r')
result_file    = open('results.fasta', "w")

my_dictionary = {} # fasta IDs are keys, value can be anything.
for line in my_id_file:
    my_dictionary[line[:-1]] = 'value'

#expected lenght of sequences in my fasta
length = 10

#pprint.pprint(my_dictionary)
for seq_record in SeqIO.parse(my_fasta_file, "fasta"):
    if seq_record.id in my_dictionary:

        ########################################################################################
        #use this line to write standard fasta record using SeqIO.write

        # SeqIO.write(seq_record, result_file, "fasta-2line") # writes fasta in 2 line format

        # SeqIO.write(seq_record, result_file, "fasta") # writes fasta with width of 60 as standard biopython format

        ### NOT POSSIBLE !!!!!!!!!!!!
        # seq_record.seq = textwrap.fill(seq_record.seq , width=10) ### --> AttributeError: 'Seq' object has no attribute 'expandtabs'

        # print(seq_record.seq, type(seq_record.seq), dir(seq_record.seq))

        sequence = str(seq_record.seq)

        result_file.write('>'+seq_record.id+'\n')

        while len(sequence) > 0:
            result_file.write(sequence[:length]+'\n')
            sequence = sequence[length:]
        result_file.write('\n')

        ### if you try to add EOL inside seq_record.seq of course you get
        ### TypeError: SeqRecord (id=seq1) has an invalid sequence.
        # SeqIO.write(seq_record, result_file, "fasta") # writes fasta with sequence width as 60   !!!! TypeError: SeqRecord (id=seq1) has an invalid sequence.

my_fasta_file.close()
my_id_file.close()
result_file.close()

```

input files :

`fasta_id_records.txt`:

```
seq1
seq2
seq3

```

`fasta_file.fasta`:

```
>seq0
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
>seq1
TTTT
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
TTTTTTTTTTTTTTT
>seq3
CCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCGGGGGGGG
GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG

```

output file:

`v` :

```
>seq1
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTT

>seq3
CCCCCCCCCC
CCCCCCCGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGCCCCC
CCCCCCCCCC
CCCCCCCCCC
CCCCCCCCCC
CCCCCCCCCG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
G

```

---
Tags: python, fasta, biopython
---


# Question
Title: How to get single histogram plot for all groups in different colours?
I am using the below script to get the histogram plot for all values together. However, I want a single plot but different colours for different groups (A, B, C)

`data file`

```
Group_name  code-1    code-2    Expression
GroupA  6  6  1.21
GroupA  6  8  0.611
GroupA  6  1  0.448
GroupA  6  3  0.3587
GroupB  8  6  0.611
GroupB  8  8  2.11
GroupB  8  1  0.60113
GroupB  8  3  0.534858
GroupC  1  6  0.448
GroupC  1  8  0.60113
GroupC  1  1  2.2
GroupC  1  3  0.622206

```

`Script`

import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams["figure.figsize"] = (12,8)

g1 = sns.histplot(df.Expression)
g1.set(
    xlim=(-2,2),
    title="Expression"
)
g1.axvline(0.0, c='black', linestyle="--")
g1.axvline(df.Expression.median(), c='red', linestyle="--")

```
</code></pre>

# Answer
> 1 votes
Just change your `g1 = sns.histplot(df.Expression)` to `g1 = sns.histplot(data=df, x="Expression", hue="Group_name")`

Here's the full code:

```
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd

plt.rcParams["figure.figsize"] = (12,8)

df = pd.DataFrame({
    "Group_name": ["GroupA", "GroupA", "GroupA", "GroupA", "GroupB", "GroupB", "GroupB", "GroupB", "GroupC", "GroupC", "GroupC", "GroupC"],
    "code-1": [6, 6, 6, 6, 8, 8, 8, 8, 1, 1, 1, 1],
    "code-2": [6, 8, 1, 3, 6, 8, 1, 3, 6, 8, 1, 3],
    "Expression": [1.21, 0.611, 0.448, 0.3587, 0.611, 2.11, 0.60113, 0.534858, 0.448, 0.60113, 2.2, 0.622206]
})

g1 = sns.histplot(data=df, x="Expression", hue="Group_name")
g1.set(
    xlim=(-2,2),
    title="Expression"
)
g1.axvline(0.0, c='black', linestyle="--")
g1.axvline(df.Expression.median(), c='red', linestyle="--")

g1.figure.savefig("hist.png")

```

'''

---
Tags: python, visualization
---


# Question
Title: Is there a tool to visualize RNA secondary structure?
I want to manually fold and pair two sequences (U4 and U6). And I want to make niceish figures.

I haven't found a single program that does that.

thanks

# Answer
The RNA viewer that is particularly good is VA RNA. This stands for "Visualisation Applet for RNA". Its a java app and run it via,

```
java -jar /Route2Progam/vaRNA/VARNAv3-93-API.jar

```

Where `Route2Progam` is the path to where the jar files is stored, all other info in the run code is literal.

Very nice pictures, point and clicky, drag structures into different shapes and easy to colour ... so for a GUI its great.

They recommend PsuedoViewer for psuedoknots, but if that doesn't mean much ... it doesn't matter.
> 1 votes

# Answer
You could try mxfold2. https://github.com/mxfold/mxfold2

```
% mxfold2 predict test.fa
>DS4440
GGAUGGAUGUCUGAGCGGUUGAAAGAGUCGGUCUUGAAAACCGAAGUAUUGAUAGGAAUACCGGGGGUUCGAAUCCCUCUCCAUCCG
(((((((........(((((..((((.....))))...)))))...................(((((.......)))))))))))). (24.8)

```
> 1 votes

---
Tags: visualization, rna
---


# Question
Title: How to remove empty files from a directory?
I have a directory of FASTA files, some of which are empty (placeholders that were output from a pipeline) and I would like to delete the empty files using a for loop on the command line.

From what I have seen so far I can use the `wc -l [file]` approach to see "if there are lines present, delete the file" but I want to make sure the syntax is right before deleting important files. I was thinking it might be easier for the command to look for the " \> " symbol found at the header in FASTA files and if it's not present delete the file.

Here is what I have so far:

```
for file in *.fasta; do "if wc -l < 1 rm file"; done

```

Or

```
for file in *.fasta; do "grep > | rm file"; done

```

# Answer
This isn't the most efficient approach, using `find` will be easier and faster as long as your `find` supports the `-empty` option, but since this is what you asked for, here is how to do it based on the number of lines:

```
for file in *.fasta; do
 [[ $(wc -l < "$file") -eq 0 ]] && echo rm -- "$file"
done

```

If you want to recursively descend into subdirectories, use:

```
shopt -s globstar
for file in  *.fasta */**.fasta; do.
  [[ $(wc -l < "$file") -eq 0 ]] && echo rm -- "$file"
done

```

If that prints the right files, remove the `echo` to actually delete them.

Note that this will also delete files that are not empty but don't have a newline character. These are not valid "text files", but strictly speaking they aren't empty, they just have no *lines* as lines are defined by the `\n` character. For example, this file would be deleted:

```
printf 'aaa' > file

```
> 2 votes

# Answer
Always Google before asking on forums. I googled "delete empty files linux" and this was the second result: https://www.baeldung.com/linux/delete-empty-files-dirs

```
# Print and delete empty files just in the current directory
$ find . -maxdepth 1 -type f -empty -print -delete
# Print and delete empty files in the current directory and all sub-directories
$ find . -type f -empty -print -delete

```
> 3 votes

# Answer
As already stated, one way is to use the GNU `find` (the default on Linux) command. This will find and print all empty files in the current directory:

```
find . -maxdepth 1 -type f -empty -print

```

Append the `-delete` action to the command above when you're ready to remove them. If you need to find (and then delete) empty files recursively, simply remove the `-maxdepth 1` option.

---

---

In a comment, you mention that these are actually "alias" files which are the results of a Nextflow pipeline. Unless the default publishDir mode has been changed, Nextflow will publish output files to the publishDir using symbolic links to the files in the working directory (i.e. `./work`). I suspect what you are trying to do is to find and remove the links to these empty files, but not the files themselves. In this case, I think what you want is:

```
find . -type l -exec test ! -s '{}' \; -print

```

Again, append the `-delete` action to the command above to actually remove the symlinks. If you instead need to delete files from Nextflow's working directory, use `nextflow clean` in conjunction with `nextflow log` to safely clean up the failed runs.
> 2 votes

# Answer
@Ram RS has identified the solution here via `find`, if that does not work then the following will work,

```
find . -maxdepth 1 -type f -empty | xargs rm

```

---

**Edit 2**

The solution is obvious ...

```
find . -type l -delete

```

The code under **Edit** is is better if the OP knows in advance *precisely* how many links the symbolic link file has, basically its always one. However, we don't know, thus `-type l` will delete every symbolic link.

---

To address @Steve's point of inadvertent deletion in case of a null return ... The code was checked prior release and this didn't occur (and also now rechecked).
> 1 votes

---
Tags: fasta, bash, grep, command-line
---


# Question
Title: Batch correction in differential expression analysis
I have sent two sets (two batches in matter of sending for sequencing) of different samples (plasma) to small RNA-seq to Qiagen company

This is how my meta data look

```
> head(meta)
           batch group    sex
30678.009 batch2     C   Male
30678.010 batch2     C Female
30941.001 batch1     C Female
30941.002 batch1     C   Male
30941.003 batch1     C   Male
30941.004 batch1     C   Male
> tail(meta)
           batch group    sex
30941.016 batch1    LT Female
30941.017 batch1    LT   Male
30941.019 batch1    LT   Male
30941.022 batch1    LT Female
30941.023 batch1    LT   Male
30941.024 batch1    LT Female
>

```

I want to extract differentially expressed RNAs between my two groups without batch effect

I am doing this

Is this correct please?

```
dds <- DESeqDataSetFromMatrix(countData = count,
                              colData = meta,
                              design= ~ batch + group)

dds <- DESeq(dds, test="LRT", reduced=~batch)

res <- results(dds)
resultsNames(dds)

[1] "Intercept"              "batch_batch2_vs_batch1" "group_LT_vs_C"
>

```

# Answer
If both batches contain replicates of both groups (LT and C) then the presented code is what you should do, exactly pointed out in the vignette for the LRT: http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#likelihood-ratio-test
> 2 votes

---
Tags: r, deseq2, differential-expression, edger, batch-effects
---


# Question
Title: How to identify a low proportion cell subpopulation in the single-cell RNA-seq data?
The cell subpopulation that I am interested in only accounts for around 1.2% of the total cells. I have previous FACS experiments that sort out the subpopulation from the samples (using markers CD166+ CD49fhi CD104- Lin-, see images). May I ask if there are any ways to identify this cell population from single-cell RNA-seq data of the tumor epithelial samples using the marker information?

# Answer
> 2 votes
Run the usual steps including clustering and visualization via something like UMAP and then color the UMAP by these four markers. That assumes that the surface marker separation you see in flow holds true on a transcriptional level which is not necessarily the case. For example, hematopoietic surface-defined populations such as the Granulocyte-Macrophage Progenitor (GMP) are well-separated in flow but on transcriptional single-cell cell level span many clusters.

Alternatively, check the tool `UCell` which can calculate a score per cell based on a combination of positive and negative markers. You can then color the UMAP by this to see whether a single cluster stands out or whether this celltype contributes to many clusters ir is not even present at all. Include some marker combinations that you know are oresent and some nonsense combination that must be absent to have positive and negative controls. That will give you a feel whether that score for your cells is reliable.

This all assumes that the cells are even in the dataset. Depending on cell size and frequency (underrepresentation) they might even have dropped out during library prep.

---
Tags: scrnaseq, single-cell, clustering, cancer
---


# Question
Title: MiXCR: only create a single export file for all clonotypes
I am using following command from MiXCR to both align, assemble and export my input files:

```
mixcr analyze amplicon \
    -s hsa \
    --starting-material rna \
    --5-end no-v-primers \
    --3-end c-primers \
    --adapters adapters-present \
    --export "-aaFeature CDR1 -nFeature CDR1 -aaFeature CDR2 \
              -nFeature CDR2 -aaFeature CDR3 -nFeature CDR3 \
              -aaFeature VGene -cGenes -count" \
    --assemble "-OassemblingFeatures=[CDR1,CDR2,CDR3]" \
    INPUT_FILE.fasta analysis_name

```

Unfortunately, this creates multiple output file (analysis\_name.clonotypes.ALL.txt, analysis\_name.clonotypes.IGH.txt, etc.), which stands in contrast to the output of the `mixcr exportAlignments` command.

How can I export a single txt file instead, while still using the `mixcr analyze` command?

# Answer
> 3 votes
This is from the MiXCR documentation:

```
Under the hood analyze amplicon is equivalent to the execution of the following
MiXCR actions options:
> mixcr align -s hs -p kAligner2 input_R1.fastq input_R2.fastq alignments.vdjca
... Building alignments
> mixcr assemble alignments.vdjca clones.clns
... Assembling clones
> mixcr exportClones --chains IGH clones.clns clones.txt
... Exporting clones to tab-delimited file

```

So when using `analyze amplicon`, you are using `mixcr exportClones` under the hood and not `exportAlignments`, hence the difference in output.

I believe what you are experiencing results from a bug, please see this issue. If this is the case, all these files should be the same according to the bug report. Fortunately the developers have already addressed the bug in version 3.0.3.

# Answer
> 3 votes
Currently in MiXCR SERIES 4 you can use `--dont-split-files`  parameter here. Please check our new documentation portal: https://docs.milaboratories.com/mixcr/reference/mixcr-export/#clonotype-tables

Also now we support presets that allow to perform analysis for a lot of standard protocols using a one line command.

# Answer
> 2 votes
After finding the correct entry in the MiXCR documentation (thanks to haci), I learned that the parameter `--receptor-type` is responsible for which sequences will be returned.

The default for this parameter when using mixcr analyze is `--receptor-type xcr`, and stands for 'all chains', thus returning such a vast output. Changing this parameter to e.g. BCR solved my issue.

---
Tags: sequence-alignment, fasta, bash
---


title: Mixcr batch processing

text:
# Question
Title: Mixcr batch processing
I have a folder with 150 pairs of fastq files from an illumina sequence run. I need to precess them with mixcr, how can I do this in Bash with a single command?

# Answer
> 4 votes
Here's how I'd approach this:

```
for READ1 in /Documents/Data/Input/*R1*
do
    READ2=${READ1/R1/R2}
    mixcr analyze amplicon --s hs \
        --starting-material dna \
        --5-end v-primers \
        --3-end j-primers \
        --adapters no-adapters \
        --receptor-type trb \
        --only-productive "$READ1" "$READ2"
done

```

Basically, loop over all the R1 FASTQ files in a directory, and calculate the R2 names as needed.

This is using BASH string manipulation, and assumes that substituting `R2` for `R1` will correctly calculate the R2 filename from R1. It also assumes that all the FASTQ files are in a single directory.

# Answer
> 2 votes
I would also suggest using GNU parallel as described in our guides :

https://docs.milaboratories.com/mixcr/guides/generic-multiplex-bcr/#one-command-solution

Briefly:

```
ls /Documents/Data/Input/*_R1* | \
    parallel -j 4 \
    'mixcr analyze
    …
    {} \
    {=s:R1:R2:=} \
{=s:.*/:/Documents/Data/result/:;s:_R.*::=}'

```

That will take all fastq pairs from /Documents/Data/Input and put the result output files into /Documents/Data/result folder.

---
Tags: ngs
---


# Question
Title: Merge data frame on either of the 2 columns in R or python
I have two data frames. One data frame with one column and second data frame with two columns. I need to merge first data frame with either column of the second data frame and returns the values. After merging the values,the order of the first column should be same as the input. Please find the sample input below

```
first_column = c("ENSG00000165588","ENSG00000213551")
df1 = data.frame(first_column)
first_column = c("ENSG00000142192", "ENSG00000140575", "ENSG00000165588", "ENSG00000165588", "ENSG00000213551", "ENSG00000213551","ENSG00000197153")
second_column = c("ENSG00000165588", "ENSG00000165588", "ENSG00000186908", "ENSG00000135446", "ENSG00000273983", "ENSG00000274267","ENSG00000213551")
df2 = data.frame(first_column,second_column)

```

Please find the output below.

```
first_column = c("ENSG00000165588","ENSG00000165588","ENSG00000165588","ENSG00000165588","ENSG00000213551","ENSG00000213551","ENSG00000213551")
second_column = c("ENSG00000142192","ENSG00000140575","ENSG00000186908","ENSG00000135446","ENSG00000273983","ENSG00000274267","ENSG00000197153")
output = data.frame(first_column,second_column)

```

# Answer
I think this question is more about subsetting than merging. You can try this solution in R:

```
first_column = c("ENSG00000165588","ENSG00000213551")
df1 = data.frame(first_column)
first_column = c("ENSG00000142192", "ENSG00000140575", "ENSG00000165588", "ENSG00000165588",
                 "ENSG00000213551", "ENSG00000213551","ENSG00000197153")
second_column = c("ENSG00000165588", "ENSG00000165588", "ENSG00000186908", "ENSG00000135446",
                  "ENSG00000273983", "ENSG00000274267","ENSG00000213551")
df2 = data.frame(first_column,second_column)

output = rbind(df2[df2$first_column %in% df1$first_column,],
               df2[df2$second_column %in% df1$first_column,])

output <- transform(output, second_column = ifelse(output$second_column%in%df1$first_column,
                                                   first_column, second_column),
                            first_column = ifelse(output$second_column%in%df1$first_column,
                                                  second_column, first_column))
output <- output[ order(row.names(output)), ]

```
> 1 votes

# Answer
**Edit** The OP in hindsight is requesting an index merge. The first solution, presented here as **Column merge solution** was rejected by the OP.

From the comments below this post the order requested by the OP was:

```
# 0 ENSG00000165588 ENSG00000142192
# 1 ENSG00000165588 ENSG00000140575
# 2 ENSG00000165588 ENSG00000186908
# 3 ENSG00000165588 ENSG00000135446
# 4 ENSG00000213551 ENSG00000273983
# 5 ENSG00000213551 ENSG00000274267
# 6 ENSG00000213551 ENSG00000197153,

```

**Column merge solution**

* The first operation is a concatenation and inner merge *****.
* The second operation are two outer merges, the second requiring a 'left-hand' merge *****.

The code is here:

```
import pandas as pd

df1 = pd.DataFrame({'A':["ENSG00000165588","ENSG00000213551"]})
df2i = pd.DataFrame({'A':["ENSG00000142192", "ENSG00000140575", "ENSG00000165588",\
                          "ENSG00000165588", "ENSG00000213551", "ENSG00000213551", "ENSG00000197153"]})
df2ii = pd.DataFrame({'A':["ENSG00000165588", "ENSG00000165588", "ENSG00000186908", \
                            "ENSG00000135446", "ENSG00000273983", "ENSG00000274267","ENSG00000213551"]})

# Column 'A'
df3in = pd.concat([df2i, df2ii])
df4inner = pd.merge(df3in,df1, on='A', how='inner')

# Column 'A+'
df3out = pd.merge(df2i, df2ii, on='A', how='outer')
df4outer = pd.merge(df3in,df1, on='A', how='outer', indicator = True).query('_merge=="left_only"')
df5 = df4outer.drop('_merge', axis=1).reset_index(drop=True).rename(columns={'A':'A+'})
df6 = df4inner.join(df5)

print (df6)

```

**Output**

```
                 A               A+
0  ENSG00000165588  ENSG00000142192
1  ENSG00000165588  ENSG00000140575
2  ENSG00000165588  ENSG00000197153
3  ENSG00000165588  ENSG00000186908
4  ENSG00000213551  ENSG00000135446
5  ENSG00000213551  ENSG00000273983
6  ENSG00000213551  ENSG00000274267

```

This approach was rejected (comments) because of the position of `ENSG00000197153` at row 2, Column A+ in relation to `ENSG00000165588`.

---

**Index merge solution**

```
import pandas as pd

def indexmerge(dfa, dfb, df1):
    df3out = dfa.reset_index().merge(dfb, how="left").set_index('index')
    df4outer = pd.merge(df3out,df1, on='A', how='outer', indicator = True).query('_merge=="left_only"')
    return df4outer

if __name__ == '__main__':

    df1 = pd.DataFrame({'A':["ENSG00000165588","ENSG00000213551"]})
    df2i = pd.DataFrame({'A':["ENSG00000142192", "ENSG00000140575", "ENSG00000165588",\
                              "ENSG00000165588", "ENSG00000213551", "ENSG00000213551", "ENSG00000197153"]})
    df2ii = pd.DataFrame({'A':["ENSG00000165588", "ENSG00000165588", "ENSG00000186908", \
                                "ENSG00000135446", "ENSG00000273983", "ENSG00000274267","ENSG00000213551"]})
    df3in = pd.concat([df2i, df2ii])
    df4inner = pd.merge(df3in,df1, on='A', how='inner')

    df4outer = indexmerge(df2i,df2ii,df1)
    df4outer2 = indexmerge(df2ii,df2i,df1)
    df4outer2.rename(columns={'A':'_A'},inplace=True)
    df6 = pd.concat([df4outer,df4outer2],axis=1)
    df6['A+'] = df6['A'].astype(str).replace('nan','') + df6['_A'].astype(str).replace('nan','')
    df7 = df6.sort_index().reset_index()
    df8 = df4inner.join(df7['A+'])
    print (df8)

```

**Output**

```
                 A               A+
0  ENSG00000165588  ENSG00000142192
1  ENSG00000165588  ENSG00000140575
2  ENSG00000165588  ENSG00000186908
3  ENSG00000165588  ENSG00000135446
4  ENSG00000213551  ENSG00000273983
5  ENSG00000213551  ENSG00000274267
6  ENSG00000213551  ENSG00000197153

```

This solution is identical to the order requested in the comments and is now also presented at the beginning of this post. The solution used here pivots around line 1 of `indexmerge()` and thus is an index merge operation. The remaining code is common `pandas` manipulation supporting this operation.

**Note** the new position of `ENSG00000197153` row 6 column A+. The index merge retains `ENSG00000197153` at the base of the column and therefore in association, i.e. same row, as `ENSG00000213551`.

---

***** technical `pandas` terminology
> 0 votes

---
Tags: r, python, dataframe, pandas, merge
---


title: TCR-seq or scRNA-seq

text:
# Question
Title: TCR-seq or scRNA-seq
I need idea, intuition, suggestion please

I have `10X scRNA-seq` of `PBMC (multiome's 3' poly-A capture)`, can I capture `TCR/BCR` from my scRNA-seq or I need `TCRseq` ? I googled and likely TCR and BCR immune profiling require 5' transcript capture and the V(D)J variable sequence is on the 5' end of the gene. So the 3' capture typically yields the C constant region of TCR and BCR transcripts. can you help me please? Thank you

# Answer
The 3' gene expression protocol will capture TCR and BCR mRNAs but this may not be very helpful to you. As you already mentioned only the 3' end will be sequenced, which are the constant regions. With this data you may potentially be able to discriminate `TRAC` from `TRBC1` and so forth.

However, in most cases researchers are interested in clonotypes (combinations of variable regions) and these you can only obtain by using the dedicated TCR or BCR kits from 10X. They contain specific sets of primers for these genes to cover the whole TCR or BCR gene body and not only the 3' end. This way you can see the different clonotypes in your population

The human TCR loci can be found for example at HUGO. The last gene at each locus is the constant and is identified by a C at the end. TRAC = T cell receptor alpha constant (I did not use the correct naming convention in the original answer, sorry)

B-cell receptors basically are membrane bound immunoglobins (IG). Therefore you will likely find the constant regions of the heavy chain in your data. The different heavy chains are what identifies the IG. For example IGHA1 is the heavy chain of an alpha IG. If the secreted and membrane bound form of the IGs is always the same per B-cell, unfortunately, is beyond my knowledge of B-cells.
> 5 votes

# Answer
According to our experience 3’ sequencing does not produce enough data to recover TCR and BCR repertoire. But we have some experience recovering repertoires from non enriched 5’ 10x RNA-Seq GEX data.

With MiXCR we have specific presets that allow analyse 10x data in one single line command:

```
mixcr analyze 10x-5gex-full-length \
 --species hsa \
 sample_R1.fastq.gz \
 sample_R2.fastq.gz \
 sample_result

```

Also feel free to check out our documentation portal:

https://docs.milaboratories.com/mixcr/reference/overview-built-in-presets/#10x
> 2 votes

---
Tags: scrnaseq, seurat, 10x-genomics
---


title: Creating a count table using R

text:
# Question
Title: Creating a count table using R
I’m trying to develop a table from a series of lists generated using the Immunarch library to process TCR sequencing results. Each list is comprised of CDR3.aa (clone) information which are character strings and their count in a particular sample. The clones are short but vary between 7 and 20 (or more). Each list has a header that identifies the sample sequenced. I have a list of 66 samples. Each sample can have several thousand clone strings. Not every clone is contained in every sample so the number of clones listed in the samples varies. Here’s an example of the structure of a single sample list.

```
$PreAg_18_2
CDR3.aa           Count
CASSYGTAYTGELFF   1623
CASSRGDSDNSPLHF   1440
CASSREKAFF        1161
CSGMGALAKNIQYF     949
CSAYTGLSYEQYF      813
CASSLSLAVNSPLHF    634
CAIRDTPGSPQHF      574
CATGQVNTEAFF       555
CASSLKGQGGSPLHF    499
CASSYSRSPQPQHF     478

```

I want to combine the results in a single table showing the clone counts with all the clone strings listed on the y-axis and the sample ID listed on the x-axis. For example:

```
                10_pep_10_1     preAg_10_2      Dec_2_18_1  …...
CASSYGTAYTGELFF    1623         234             0
CASSRGDSDNSPLHF    1440         522             28
CASSREKAFF         1161         445             50
CSGMGALAKNIQYF      949         24              0
CASSYSRSPQPQHF      478         0               398
.
.

```

Currently, I'm working with R because that's what Immunarch is in, but I wonder if a) There's a better way to do this in Python, or b) If it's doable in r, then how do I generate the matrix matches. The main problem is that not all the samples have the same number of clones, so I first have to separate each sample data before I can convert it to a data frame, and then how can I match the actual count with the actual clone. I can't just merge the data frames. Any suggestion is greatly appreciated.

Here's what i have for code right now. Code: library(immunarch) library(stringr) library(plyr)

```
immdata = repLoad("/mnt/data/Development/Analysis_Scripts/input_files/")

all <- immdata$data

# Get list headers (names) and convert to df
sample.id <- names(all)
sample.id <- data.frame(sample.id)

# Get list of clones and filter for unique clones per list.
for (i in 1:length(all)){
    all[[i]]$Sample.ID<-names(all)[i]
    all[[i]]<-all[[i]][,c("CDR3.aa", "Clones")]
}

# all is the variable that contains the list of samples and clones.
all.split <- split(all, )

# make vector of all clones
all.clones <- unlist(all, use.names=FALSE)

# Removes clone repeats
all.clones.u <- unique(all.clones)

# convert list of clones to data frame
all.clones.u <- data.frame(all.clones.u)

```

At this point I have both the list of sample ids and the list of clones. My problem is now matching each count with respective sample and clone in the table. Any suggestions?

# Answer
> 4 votes
You could give a try with dplyr's `full_join()` and purrr's `reduce()`. `full_join()` will make sure that you won't lose any data and `reduce()` will merge tables in your list from left to right, left being the first table of the list in the first iteration but would become merged *table1+table2* in the second iteration and so on. This left to right "reduction" will also be helpful when renaming the columns of the resulting table.

Without your workspace I have attempted to replicate your data: There is a list of 3 data frames, each with 2 columns:

```
library(dplyr)
library(purrr)

df1 <- data.frame(CDR3.aa = c("A", "B", "C"),
                  Count = c(10,20,100))

df2 <- data.frame(CDR3.aa = c("A", "B", "D"),
                  Count = c(100,200,300))

df3 <- data.frame(CDR3.aa = c("D", "E", "F"),
                  Count = c(10,20,100))

my_list <- list(df1 = df1, df2 = df2, df3 = df3)

> my_list
$df1
  CDR3.aa Count
1       A    10
2       B    20
3       C   100

$df2
  CDR3.aa Count
1       A   100
2       B   200
3       D   300

$df3
  CDR3.aa Count
1       D    10
2       E    20
3       F   100

# "by" argument specifies which column to use when merging
my_table <- reduce(my_list, full_join, by = "CDR3.aa")

names(my_table) <- c("CDR3.aa", names(my_list))

> my_table
  CDR3.aa df1 df2 df3
1       A  10 100  NA
2       B  20 200  NA
3       C 100  NA  NA
4       D  NA 300  10
5       E  NA  NA  20
6       F  NA  NA 100

```

Clonotypes that are not present in table will be "NA".

# Answer
> 1 votes
That can be easily done with MiXCR software. If you first analyze the raw data using mixcr, which can be easily done with `mixcr analyze` command depending on your data structure:

e.g.

```
mixcr analyze generic-tcr-amplicon \
--species mmu \
--rna \
--rigid-left-alignment-boundary \
--floating-right-alignment-boundary C \
  input_R1.fastq.gz \
  input_R2.fastq.gz \
  result

```

More on that here:

https://docs.milaboratories.com/mixcr/reference/mixcr-analyze/

After that you will have .clns file for every sample.

Then, to merge all samples you can use `mixcr exportClonesOverlap` and all .clns files as an input, e.g.:

```
mixcr exportClonesOverlap \
  --criteria "CDR3|AA|V|J" \
  clonesets/*.clns
  overlapTable.tsv

```

`--criteria` allows you to specify the criteria required for clonotypes to be treated as identical. `"CDR3|AA|V|J"` means that two clonotypes have to share cdr3 amino acid sequence, v gene and j gene.

overlapTable.tsv is exactly what you need.

There are other parameters you can tune (like export only productive clonotypes etc.), more on that here:

https://docs.milaboratories.com/mixcr/reference/mixcr-exportClonesOverlap/#example

---
Tags: r
---