jlsteenwyk / clipkit Goto Github PK

Hello.

I am having the following issue, and below I post how I fixed it. Thank you for the great tool.

The issue occurred on Python 3.8 and clipkit version 1.1.3, BioPython 1.79 and Numpy 1.20.3, and using the following command line:

clipkit clean.aln -m kpic-smart-gap -l -o filtered.aln -if fasta -of fasta

Which in the past has worked well, but today my output from a FASTA alignment (requesting a FASTA output alignment) is returning numerical data:

>MN908947.3
676888886866667666888666686878787768786686776876687688678766
686667667868668866866686688668786788766676666768666867886686
876676876886678886786677887666768686676666868677888678667778
787666766677867687767676688786687788866676766666666678666686
678887668788886667788676766787687866787768887767668667877677

I have tracked the issue down to this line of code:

Essentially, that is returning an int because the code is indexing a bytes object (https://docs.python.org/3/library/stdtypes.html#bytes-objects). I am unsure if there has been some change to BioPython or Numpy here that might be causing the issue.

Changing the line to:

            keepD[entry.id][alignment_position] = entry.seq._data.decode()[alignment_position]

Makes thing work, but at a significant cost to speed.

So, I tried a different approach that works at a good speed but requires the following changes:

1. Instantiate np.zeros rather than np.empty, with dtype=bytes to store the state of locations to keep and discard:

So, it looks like this:

keepD[entry.id] = np.zeros([alignment_length], dtype=bytes)

It has to be zeros so that we get an array of b'', which will crucial below when joining things (if you use np.empty, you get random bytes that don't decode back to strings and can't be joined).

2. These two lines need to change to:

trimD[entry.id][alignment_position] = entry.seq._data[alignment_position:alignment_position+1]

This ensures you get a bytes object of length 1 back rather than an int.

3. Finally, this bit of code needs to decode the bytes arrays to strings before joining:

    for k, v in keepD.items():
        keepD[k] = "".join(np.char.decode(v))
    for k, v in trimD.items():
        trimD[k] = "".join(np.char.decode(v))

Thanks again.

Hello,

I used pip install clipkit to install and am using clipkit v0.1.

I was initially getting an error that the input file cannot be read, so I altered all files to chmod 777 just to be sure that wasn't an issue.

(ref_creation) swyka@taxon-pipeline:~/references/2020-09-10_updates$ ls -al
total 276
drwxrwxr-x  5 swyka swyka   4096 Sep 14 13:49 .
drwxrwxr-x  3 swyka swyka   4096 Sep 10 11:27 ..
-rwxrwxrwx  1 swyka swyka 129481 Sep 14 13:49 EOG09261EM7.fa
-rwxrwxrwx  1 swyka swyka 129481 Sep 14 13:47 EOG09261EM7.fasta

However, when running I am still getting the same error.

(ref_creation) swyka@taxon-pipeline:~/references/2020-09-10_updates$ clipkit EOG09261EM7.fasta -o EOG09261EM7_trimmed.fasta
Traceback (most recent call last):
  File "/home/swyka/miniconda3/envs/ref_creation/bin/clipkit", line 8, in <module>
    sys.exit(main())
  File "/home/swyka/miniconda3/envs/ref_creation/lib/python3.6/site-packages/clipkit/clipkit.py", line 111, in main
    execute(**process_args(args))
  File "/home/swyka/miniconda3/envs/ref_creation/lib/python3.6/site-packages/clipkit/clipkit.py", line 59, in execute
    input_file, file_format=input_file_format
  File "/home/swyka/miniconda3/envs/ref_creation/lib/python3.6/site-packages/clipkit/files.py", line 49, in get_alignment_and_format
    raise Exception("Input file could not be read")
Exception: Input file could not be read

Although, when I try to force the file format -if fasta I get a slightly different error.

(ref_creation) swyka@taxon-pipeline:~/references/2020-09-10_updates$ clipkit EOG09261EM7.fasta -o EOG09261EM7_trimmed.fasta -if fasta
Traceback (most recent call last):
  File "/home/swyka/miniconda3/envs/ref_creation/bin/clipkit", line 8, in <module>
    sys.exit(main())
  File "/home/swyka/miniconda3/envs/ref_creation/lib/python3.6/site-packages/clipkit/clipkit.py", line 111, in main
    execute(**process_args(args))
  File "/home/swyka/miniconda3/envs/ref_creation/lib/python3.6/site-packages/clipkit/clipkit.py", line 59, in execute
    input_file, file_format=input_file_format
  File "/home/swyka/miniconda3/envs/ref_creation/lib/python3.6/site-packages/clipkit/files.py", line 34, in get_alignment_and_format
    alignment = AlignIO.read(open(inFile), file_format.value)
AttributeError: 'str' object has no attribute 'value'

This is the head of my fasta file, for reference.

(ref_creation) swyka@taxon-pipeline:~/references/2020-09-10_updates$ head EOG09261EM7.fasta 
>XP_001588528.1|GCF_000146945.2|EOG09261EM7
MYSRYVPPSKKKVTVGDQLSQAPLSPLKSSSPPPPPSPAPAIRPDASSSYARYIPPSKSNVKPDAQLNAPKLGLESPSPASKRKREDALEPGPEPVLKKAKKDKKEKSVKVHASATLDEPHADNASSEEAQEVTKKDKKQKKKKSSDDSTSSEETTENPEDIDDKRHKKVLEKREKSIKKAERRARKAAEEGRDAEDAQEPEEPVEIHNLVPLPQPEPIPELPPPSLESTLPSWLASPILVSPTTT
AQFSEVGVEAEAATVLRSKGFNEAFAVQAAVLPLLLHGTRQKPGDVLVSAATGSGKTLSYVLPMTQDISRNIVTRLRGLIVMPTRELVSQAREVCEVCSSAFSAGSRKRVKIGTAVGNEAFKVEQANLMENTYKYDPARYHEQERRKNLRWESSDAGTDDEGEPLLDDEAISPLPDHIIEPVSKVDILICTPGRLVEHLKFTPGFTLEYVKWLVIDEADKLLDQSFQQWLNVVMSSLATGQNSFPN
NRDRVRKIVLSATMTRDIGQLTSLKLYRPKLVVLAGSSAGDDGKSSHAHILPPGLVEFAVKVDDENLKPLYLMEILKGNDMIDDSKIKSDSDTDSDTDSDSDSDDSSSDSSSDKDSSEDSSSSDDSSISGSDSESKPDKVSPKSKPKFKTNPLPISHEPHGVLIFTKSNESAIRLGRLISLIHPSYTEIIGTLTSTTRSSERKASLASFSRGKLQILVASDLVSRGLDLPDLAHVINYDVPTSITN
YVHRVGRTARAGRQGHAWTLVGNSEARWFFNEVAKSEEIRRRESAKVKRIVVDARKFGEYKKETYEEALEELGQEAMAIRSKK
>TGO71980.1|GCA_004786205.1|EOG09261EM7
MYSRYVPPSKKKVATQEKSIEAPPSPPKSSSPPPPLAPPATRPDASSTYARYIPPSKSKSKPNTQLDAPITAPESPSPASKRKHEDVVQPTPEPISKKARKEKKEKSVTALMSSVSEESHKDTISSEEAQEVTKKDKKSKKSKSKVETIPSRDTAGSNEDIDDTRHKRVLEKREKSIKKAERRARKAAEEGVAVENAQEPEEPAEIHDLVPLPQPEPIPELPPPSLESTLPSWLASPILVSPTTTA
EFSDLGVEVEAANVLRSKSFNEAFAVQAAVLPLLLPGSQQRPGDVLVSAATGSGKTLSYVLPMTQDISRNTVTRLRGLIVMPTRELVSQAREVCEVCSSAFSVGSRKRVKIGTAVGNEAFKVEQANLMENTHRYDPILYHEQEQRKNSRWESSDAGTDDEGEPFFDDEIVSPLPDHVIEPVSKVDILICTPGRLVEHLKSTPGFTLQHVKWLVIDEADKLLDQSFQQWLDIVMNSLAAGQKSLPSN
KDRVRKIVLSATMTRDIGQLTSLKLYRPKLVVLEGSSAGDDGKGSQAHVLPSGLAEFAVKVDDENLKPLYLIEILKGNNMIDESEIKSDTDTDSDSDSSSDSSSDEDSSEDSSSSDDSSESDKDSDSKPDSRPSSKSKSKIKTNPLPINHEPHGVLIFTKSNESAIRLGRLLSLINDSYTNIIGTLTSTTRSSERRASIASFSRGKLQILVASDLVSRGLDLPDLAHVINYDVPTSITNYVHRVGR
TARAGRQGSAWTLVGNTEARWFFNEVAKSEEIRRRDGEKVKRVVVDARKFGEYKKETYEDALEELGQEAMAIRKKKSS
>EPQ63759.1|GCA_000418435.1|EOG09261EM7
MNSYSYARYVPPPKSKISESIRDSTSSQPISSSKSVITCDTYHDASASYARYLPSKPISNGANSAPKFDRNSSPSSSLAKIGTRSLKPEEGPLPNIKSMNAPLSRKSPTQQGRSPSDQISSKSKRYVKTFKSDTNAQVIMDNQAYSEPSIDINKINRNKNNSNPDQRTKFNQEGLKSPSKSGSPPSGTYHQNKSEPSATIEQGKISKLDLQNVGSINVKKRARDELDSGDEKSTKNHNKILKKREK
SLKKIKAKEDASKGLPLMEPSTKVSVELHDLVPLPQAELLPPATLVSIADSYPPWMANPTHVTTQKKVSFKELGLENHTIKVLQEKGFNEAFAIQAAVLPLLISNTERERGDIVISAATGSGKTLAYTLPMIKDISYHKITKLRGLIILPTRELVFQVKEMAETCVTAFYHNKKKRVKIGTAVGNESLSVEQSHLIDYDQKYDPEEYSKRLNRIDLAWKSAEDVDDGFHSIFEESCTGKFPDHISY
QKVNVDILICTPGRLVEHMKCTSGFSIKDVSWLVIDEADKLLDQSFQQWLPLVMAEVESKHTEPLVKRRRVQKLILSATMTRDLGELAQLKLYRPKLITLDSEMEIDGDSQTHNLPCTLWESGVKVEEDGIKPLYLMEILKRESVMSIADLEDSYPSDTDSLDSSLDSETSDSSTTVSKERSTSKVLSLSPDSNSNNPQGVLIFTNSNVTAVRLGRIIPLLCPNLATRIGVLTSSLPRSSRQQCIR
SFCKGIISVIVASDLVSRGLDLPNLAHVINYDIPTSVVSYIHRVGRTARAGKEGRAWTLFTATEGRWFWNEIGRSAHIHRPKGKIVRTNISKETFSEEQRRTYGLALETLEKEATGVEICKWT

Hello,

I am interested to use your software to clean some alignment. I wanted to trim colunm with more than 90% of gaps/missing data and by keeping codon.

I tried the following command

clipkit infile -m gappy --gaps 0.9 --codon --output outfile

I ahve this error

IndexError: arrays used as indices must be of integer (or boolean) type

When I used only --gaps 0.9 the gap thresold is 0.4138. I am only able to have a threshold of 0.9 by use only the argument -m gappy. (-m gappy plus --codon doesn't work).

Is it a bug known?

I am using the last version 2.2.4 installed by conda.

thanks for the help,

regards

Nicolas

I wonder if there is an option to preserve the sequence description for FASTA alignments? We store important info about the sample in the description but they are stripped by ClipKIT. I could write code to re-attach the descriptions from the original alignments but it would be ideal if ClipKIT could leave the descriptions untouched.

Thanks for the very useful software,

Edgardo

Bug Summary

Not a bug actually. But I found the trimming behavior is different from the previous version.

In the previous version, using gappy or smart_gap mode trims sites that the gappyness above a threshold. But in the latest version it trims sites equal or above the threshold.
Here are the codes I found before implementing the numpy acceleration:

However, in the latest version:

Just want to confirm you're expecting this behavioral changes, since it will produce a slightly different results from the previous version.

Steps to Reproduce

from Bio.Align import MultipleSeqAlignment
from clipkit.msa import MSA
from clipkit.settings import DEFAULT_AA_GAP_CHARS

records = ["MG--A", "M-GTC"]
pep_msa = MultipleSeqAlignment([SeqRecord(Seq(rec)) for rec in records])
gaps = 0.5

clipkit_pep_msa = MSA.from_bio_msa(pep_msa, gap_chars=DEFAULT_AA_GAP_CHARS)
clipkit_pep_msa.trim(mode=TrimmingMode.gappy, gap_threshold=0.5)
print(clipkit_pep_msa.to_bio_msa())

The previous codes do not trim any of the sites and will get:

Alignment with 2 rows and 5 columns
MG--A <unknown description>
M-GTC <unknown description>

The current codes will trim the sites that gappyness == 0.5 and get:

Alignment with 2 rows and 2 columns
MA <unknown description>
MC <unknown description>

Technical Details

• Python version: 3.10.12
• ClipKIT version: 2.1.1

Hello teacher,
I encountered a problem during installation.
I tried Conda or the source code but encountered the following error. Can you help me see why
my python was 3.8

Traceback (most recent call last):
File "setup.py", line 4, in
from clipkit.version import version
File "/home/zhangwenda/biosoft/tree/ClipKIT/clipkit/init.py", line 1, in
from .api import clipkit
File "/home/zhangwenda/biosoft/tree/ClipKIT/clipkit/api.py", line 4, in
from .clipkit import run
File "/home/zhangwenda/biosoft/tree/ClipKIT/clipkit/clipkit.py", line 9, in
from .args_processing import process_args
File "/home/zhangwenda/biosoft/tree/ClipKIT/clipkit/args_processing.py", line 5, in
from .helpers import SeqType
File "/home/zhangwenda/biosoft/tree/ClipKIT/clipkit/helpers.py", line 6, in
from .msa import MSA
File "/home/zhangwenda/biosoft/tree/ClipKIT/clipkit/msa.py", line 16, in
class MSA:
File "/home/zhangwenda/biosoft/tree/ClipKIT/clipkit/msa.py", line 96, in MSA
def is_any_entry_sequence_only_gaps(self) -> tuple[bool, Union[str, None]]:
TypeError: 'type' object is not subscriptable

Dear author,

i use the command below to extract parsimony sites,

clipkit -l -m kpi-gappy i.fa -o o.fa

but get much more sites than the really parsimony sites,

i wonder does clipkit calculate gap(N/-) as variable sites when counting parsimony sites;

I'd like to be able to keep track of CDS & protein sequences throughout the trimming process, and I don't want to trim CDS sequences in a way that would introduce frameshifts (i.e., trimming at the codon level). Gblocks can do this, but I'd like to see if this method can do a better job than Gblocks. Essentially, this involves translating the CDS sequence to protein sequence, trim the protein sequence (keeping track of which positions get trimmed), and then reverse-translating back to CDS sequence.

Alternatively, is there a way to obtain the trimmed positions from an amino acid alignment? The -c option is close, but becomes ambiguous in gappy alignments. Better would be a list of aa positions that were trimmed. I can then take that list and reconcile the untrimmed CDS with the trimmed protein sequence alignment. Perhaps this function is already possible and I am just not reading the documentation thoroughly?

Thanks!

Some of the software pinned to a version, say click or tqdm negatively interacts with many other tools that may need newer features developed for those packages.

It is exceedingly unlikely that click or tqdm will introduce API changes that would break your code, thus I would recommend to only pin versions that have a good reason to be pinned.

For all others, you should list a minimal version number (and perhaps a maximal version number if that later turns out to be necessary).

Hello,

I wanted to know the differences in terms of output MSA and subsequent HMMs that are to be constructed, when one uses smart-gap, gappy and kpic-gappy parameters to trim.

Thanks
Abhijit

In help message,
--modes
should be
--mode

clipkit --version
clipkit 1.1.3

clipkit --help

optional arguments:

  -m, --modes <smart-gap,                     trimming mode 

Thanks for this really useful library!
I'm trying to find the code you have used to compare ClipKIT with other methods (the graphs shown here).

if you have written any wrappers around trimAI, Gblocks (whose source code I cant seem to find) and BMGE, could you please share them too?

Hi Jacob & team,

With interest I'm exploring the alignments processed by ClipKIT, including how each site is classified. For the latter, I appreciate the option to generate a log file with the option -log. However, it appears that using this option in combination with specifying an output file for the trimmed alignment itself (with the -o option) does prohibit the latter, since in this case the results are written to output file 'g' (and the classifications to 'g.log'). Possibly I'm submitting the wrong command (e.g. clipkit [alignment file] -m kpi-gappy -g 0.7 -o [trimmed alignment file] -log), in that case I'd like to figure that out.

Best regards,
Jolien

Hi - I'm suggesting to use some numpy implementation instead of for loop for trimming.

Below is an example codes that I tried. I slightly modified the original function ClipKIT used for the trimming process but the concept should be the same. Note that I'm using Biopython MultipleSeqAlignment objects as input. I noticed you have already convert them to numpy array so I guess you can simply skip the conversion step. Also the trimD definition might be different from the original one since I just use it to store the index should be trimmed.

def original_clipkit(msa: MultipleSeqAlignment, gaps: int = 0.9):    
    keepD = msa[:, :0]
    trimD = []
    for idx in range(msa.get_alignment_length()):
        _, gappyness = report_column_features(msa, idx, get_gap_chars(SeqType.aa))
        if gappyness <= gaps:
            keepD += msa[:, idx : idx + 1]
        else:
            trimD.append(idx)

    return keepD, trimD

I tested it with jupyter %%timeit function and here is the results:
67.1 ms ± 1.67 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

And here is the rewrite with numpy methods.

def numpy_clipkit(msa: MultipleSeqAlignment, gaps: int = 0.9):
    infoList = [{"id": rec.id, "name": rec.name, "description": rec.description} for rec in msa]
    msa_array = np.array([list(rec) for rec in msa])
    gappyness = (msa_array == "-").mean(axis=0)
    trimD = np.where(gappyness > gaps)[0]
    keepD = np.delete(msa_array, trimD, axis=1)

    return MultipleSeqAlignment([SeqRecord(Seq("".join(rec)), **info) for rec, info in zip(keepD.tolist(), infoList)]), trimD.tolist()

And there is a significant speed up:
616 µs ± 2.73 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

Since I'm only using the gappy mode to do the trimming, I'm not sure if this implementation is applicable to other trimming mode.

Thanks for this wonderful tool! I wanted to clarify some information in the help file about the defaults for the --gap_characters argument.

Should the default for nucleotides include N as a gap character as the information under Sequence type suggests?

  Sequence type
      Specifies the type of sequences in the input file. Valid options
      include aa or nt for amino acids and nucleotides. This argument
      is case insensitive. This matters for what characters are
      considered gaps. For amino acids, -, ?, *, and X are considered
      gaps. For nucleotide sequences, the same characters are
      considered gaps as well as N.

I just wanted to check since the defaults under --gap_characters reflect the opposite I believe.

  -gc, --gap_characters <string_of_gap_chars> specifies gap characters used in input file
                                              (default for aa: "-?*XxNn"
                                               default for nt: "-?*Xx")

Thanks very much!

Technical Details

• Python version: Python 3.10.12
• ClipKIT version: clipkit 2.2.2

Hello,

while trying to install clipkit 2.11 using conda, I get the following error.

Could not solve for environment specs
The following packages are incompatible
├─ biopython 1.81** is installable with the potential options
│ ├─ biopython 1.81 would require
│ │ └─ python >=3.10,<3.11.0a0 , which can be installed;
│ ├─ biopython 1.81 would require
│ │ └─ python >=3.11,<3.12.0a0 , which can be installed;
│ ├─ biopython 1.81 would require
│ │ └─ python >=3.12,<3.13.0a0 , which can be installed;
│ ├─ biopython 1.81 would require
│ │ └─ python >=3.8,<3.9.0a0 , which can be installed;
│ └─ biopython 1.81 would require
│ └─ python >=3.9,<3.10.0a0 , which can be installed;
└─ pin-1 is not installable because it requires
└─ python 2.7.* , which conflicts with any installable versions previously reported.

Pins seem to be involved in the conflict. Currently pinned specs:
- python 2.7.* (labeled as 'pin-1')

Do I have to install both python 2.7 and python 3.12 in the same environment before proceeding?

The following command only installs version 1.4
conda install -c bioconda clipkit

Is there a lot of difference between version 2.1.1 and 1.4? If not and all else fails, I might just use version 1.4

Thanks
Abhijit

Dear Clipkit Developer

I am interested using Clipkit before running IQTree. But I also need to generate Consensus sequences for primer designing. Should I used Clipkit trimmed alignment for Consensus Sequence generation or not?

Also, you autodetect, detects my Nucleotide sequence as protein sequence:

Determining smart-gap threshold...

-------------
| Arguments |
-------------
Input file: ../05_postOutlierDetection/halign3.aln (format: fasta)
Output file: ../05_postOutlierDetection/halign3.aln.clipkit (format: fasta)
Sequence type: Protein
Gaps threshold: 0.9986
Trimming mode: kpic-smart-gap
Create complementary output: False
Create log file: True


------------------------
| Processing Alignment |
------------------------