text:
Compare VCF files
I have two strains, X and Y, and I would like to compare the variants found in both. I have my own parent strain, Z, that was sequenced and assembled to be used as the reference genome.
I first aligned X and Y to Z using minimap2 and then I used FreeBayes as the variant caller. Once I had the two separate VCF files I looked online to see what was the best way to filter out the unique and common variants in X and Y.
I found the bcftools isec tool to be my best bet so I went ahead and ran this line of code:
./bcftools isec -p out/path/ X.vcf.gz Y.vcf.gz
Now I have three new VCF files:

0000.vcf (variants unique to X)
0001.vcf (variants unique to Y)
0002.vcf (variants in both X and Y)

However, it's still in VCF format. I need a way to parse through the VCF file to list out the SNPs with high quality with their position on the genome. I tried using bcftools stats but it just lists out how many SNPs there are and not their position.
TLDR: Can someone help me find a way to list out the locations of SNPs from a VCF file in human readable format?

Can you explain what you need a bit more clearly? Try showing us a few lines of the VCF file and what you would like to convert them to, for example. Since the VCF already gives locations of variants in a human readable format, your question is a bit unclear. Do you want SNPs only, without any larger variants? Do you just want the first, second, fourth and fifth fields of the VCF? Do you want to convert it to hgvs notation? Something else?

If I understood you correctly you are searching for a way to represent the genomic locations in your vcf file. However, the vcf file is already a human readable file (vcf specification).
The first lines are the header, but afterwards you can find your information in a tab delimited fashion:
#CHROM POS ID REF ALT QUAL FILTER INFO

So the information on genomic positions is directly readable from the vcf file (first 2 columns)


text:
change the order in ggplot
I'm trying to make a plot that looks like this:

This is a rudimentary plot I did in excel. It's a plot of reads by length. Very simple nothing crazy to it. The important thing to note is that it's plotting the data in increasing order of length.
When I try to use ggplot to make this same plot in R it does something like this:

Where it orients the read names in alphabetical order not in ascending order of the lengths. The dataframe/table I have read into R is already listed in ascending length order. For example:
1   c6dd66b4-dab5-484b-92b6-c4ecdcc155b1   1570
2   ecf15d4f-15b1-4636-a97d-e0928281e311   1739
3   461f05c7-0fc6-4564-95df-6188966fe616   1790
4   e989c32c-491e-4af1-a37d-4d4f5c87a3d9   1797
5   ce94c2b3-442b-4b8c-bd59-6d7c23d5fe3d   1824
6   ed15ea16-a28d-42ca-9a07-2818d90ddd4b   1832
7   481550cd-6f87-44c7-915f-f0f3d638ffe4   1853

As you can see the lengths increase as you go down and the read names are in no particular order. Can someone help me read in the correct order to make this plot actually work in ggplot?

the technical term you want is called "ranking". The question (and the answer) will be much clearer once that is used. rank is one of the possible R methods that will do this.

You need to enforce the order you need on the read-name column. Let's say df is your data frame with read.name and length as the two columns.
This is a fairly easy operation.
df <- df[order(df$length),]
df$read.name <- factor(df$read.name, levels = df$read.name, ordered = TRUE)
ggplot(df, ...) #Your current ggplot code goes here

Hmm this gave me an error factor level [99] is duplicated but I guess I should just google that too, right? No offense but there's no need to be rude with your response. I did Google this question extensively before making this post and the order and factor methods did pop up and I did try them but didn't get it to work which is why I posted this....

My apologies @rimo I'm a bit slow on this. The solution to your data is that it needs ranking before its plotted, for example df$rank <-  rank(df$length) then if you place the df$rank into ggplot2 that should work. RamRS is using order instead.  If it was Python I'd code it for you (coding is easier than writing). It is fairly easy with R however.

Yes, you should google that too. That's how we learn. Google error messages and try stuff until you either understand the whole picture or run into something a lot beyond your understanding. The error says your read names are duplicated. In that case, I'd either make them unique by adding a new component, use rank as M__ suggests or plot length vs 1:nrow(df) in the x axis so sorting by x-axis doesn't do a thing, as your plot doesn't really need read names on an axis.

source: Stack Exchange

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

sort: votes


text:
Combine separate fasta files from a list.txt
I have about 10 fasta files that I need to pull from a giant folder of other fasta files (has about 100 files). I have a list of the files I need in a .txt file separated by line and I want to pull out the 10 fasta files listed and combine them into one big file.
I tried this command at first:
xargs < names.txt cat > names.fasta
But it gave me this error:
cat: 'read_1.fasta'$'\r': No such file or directory
For all the files except the last read in the list which it managed to add to the names.fasta .
I'm not sure what the error message means given I can find the files in my folder so they should be there. I can do this manually but I know it'll take a while to do.
EDIT: I changed the list to be separated by space and it actually pulled out all the files I needed! However, is there a way to account for reading in the names by newline instead of space??

Please always mention what operating system you are using. Your error suggests you are using Windows and have Windows line endings (\r\n) and that is essential information when trying to manipulate files and parse text.

This should work
xargs < names.txt cat > names.fasta

Or else try,
cat names.txt | xargs cat > newdir/names.fasta

What I think is within names.txt there's a typo or a missing file and that's what it's complaining about.
Perhaps try using
ls > tmp
cat names.txt | sort > tmp2
diff tmp tmp2

If there are any differences that should pick it up.

Just to address an edge case, run dos2unix names.txt first and then a modified version of @M__'s command
cat names.txt | xargs -I v_file cat v_file >> names.fasta

Based on your edit, it looks like your working platform has different line endings compared to the platform on which names.txt was last edited/saved.

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


text:
Phylogenetic Separation
I am having a bothersome issue with my current approach in generating a phylogenetic tree of yeast samples, where samples from two groups (in-house and public NCBI data) form separate clades, despite representing the same biological sample.
I have about 50 strains going into this tree that come from two different groups:

sequenced data from our in-house strains
sequenced data from NCBI (of samples where we directly obtained our in-house strains)

Simply put, sample_1 from group_1 is directly sourced from sample_1 from group_2 so in theory these two samples should sit right next to each other in a phylogenetic tree. However, my final tree shows two different clades (each with their own branching) which consist of group_1 for one clade and group_2 for the other - so they are completely separate in the tree.
I tried to troubleshoot this by:

making sure all the samples have the same read lengths (and trimming adapters)
downsample each sample to the same coverage where I can
normalizing the multi-sample final VCF file (bcftools norm)
filtering the VCF file with basic filters (bcftools filter -i 'QUAL > 30 && FORMAT/DP > 10 && MQM > 40 && (QUAL / FORMAT/DP) > 2.0')
filtering the allele frequency (AF) to be > 0.25

Doing this helped put a handful of samples next to their partner from the other group in a new tree I've made but still not all of them are together the way they should be. I assume the biggest issue comes from the difference in sequencing runs / how I calculate the distance so how would you troubleshoot this error further? Has anyone experienced an issue like this?
This post is already so long but in case it helps, this is my current approach:

pre-process the sequences as mentioned above (standardize read length and coverage, trim adapters)
align sequences to reference genome using BWA MEM
mark duplicates with picard, add read groups with picard, mark duplicates again
filter the final alignment by quality (samtools view -q 50)
call variants using freebayes and a list of every bam file to call simultaneously
use gatks VariantsToTable on the vcf file to get the GT info
use a custom R script to calculate pairwise distances
use a custom R script to plot the tree

I have some suggestions. I'm assuming of course that you are using the same reference genome across samples.

Go find the positions that lead to this phylogenetic placement. You should be able to do this by crosstabulating genotypes in your VCF file by the artefactual groups. Do they have anything in common? For example, if you have a large number of mitochondrial variants but one of these groups of samples has had mitochondrial reads filtered out somehow, you might see this kind of thing.
Directly inspect a subset of alignments from the two groups using IGV. Are there any weird artifacts? Visualizing coverage, do you see anything systematic?
Are these 50 NCBI genomes from a single study? Can you just go look at methods from a paper and see if they did anything weird?


text:
Filtering file with AWK and writing output to new file
Hi I'm trying to filter a column in a PAF file with the AWK command and write it to a new file:
awk '$12>=20 {print $0}' /PATH/FILE.paf > /OUT/PATH/FITLERED_FILE.paf
This command makes sure the 12th column has a value of at least 20 and then prints the entire contents of the file - with the filtered column - to a new file.
This works for one file but I have 90+ I need to filter through. I tried using a for loop:
for i in /PATH/TO/FILES/*_mapped.paf; do awk '$12>=20 {print $0 > "awk_{}.out"}' $i; done
Which didn't work. And then I tried using a pipe command:
cat /Users/rimo/Desktop/Pat_Assemblies/flye_mapped.txt | parallel "awk '{$12>=20 print $0}' /Users/rimo/Desktop/Pat_Assemblies/Assemblies_Flye/{} > awk_{}"
Which also didn't work. I think the problem I am having is format-based but I would appreciate it if someone could steer me in the right direction!

What operating system are you using? Do you have access to GNU awk (gawk on some systems, just awk on most Linuxes except recent Ubuntu and maybe some others).

I'm using a mac OS terminal

OK. Do you have gawk? No need to install it if you don't, it would just make some things easier if you already have it.

I don't have gawk but I could install it if it's easier! What would you suggest as an alternative?

You don't really need gawk here. I was mistakenly thinking you wanted to edit the original file which is made easier by gawk (you can use gawk -i inplace) and gawk won't choke if you have too many open files, but all you need to avoid that is to close the file handle after writing (see the last command in my answer), so you really don't need to worry. Just use the default awk.

The {} isn't in any way special to the shell or awk. It is used as a replacement by certain other programs, off the top of my head I can only think of find and xargs where it is used by convention, but it isn't a shell thing or an awk thing. Here, you are using it inside the awk command so it just breaks the awk's format because { } introduces a command block and this one is empty.
Your file name is $i, so if you want to make foo_mapped.paf into awk_foo_mapped.paf.out, you need awk_$i.out. Something like this:
for file in /PATH/TO/FILES/*_mapped.paf; do
    fileName="${file##*/}"
    awk '$12>=20' "$file" > awk_"$fileName".out
done

Note how I removed the print $0. This is because the default action for awk when something evaluates to true is to print the current line, so awk '$12>=20 {print $0 } is equivalent to awk '$12>=20.
You can also use your original approach of having awk deal with the file name by using awk's internal FILENAME variable. Like this:
awk '$12>=20{ print > "awk_"FILENAME".out" }' /PATH/TO/FILES/*_mapped.paf

Depending on the number of files, you might need to close the file handle after writing. If you get an error about too many open files, try:
mawk '$12>=20{print > "awk_"FILENAME".out"; close("awk_"FILENAME".out")}' /PATH/TO/FILES/*_mapped.paf


text:
hierarchical clustering of kmers and their counts
I have a list of kmers and their frequencies broken down like this:
TGCCAA  2
TGGACA  2
TGGGAA  2
TGGTAA  5
TGTCAA  5
TGTGAA  2
TTCCAA  13
TTGAAA  6
TTGCAA  4
TTTAAA  6

And I would like to use clustering to look at the most frequent and the most unique instances of kmers in a sample. I started off following the tutorial found here (https://www.w3schools.com/python/python_ml_hierarchial_clustering.asp#:~:text=Hierarchical%20clustering%20is%20an%20unsupervised,need%20a%20%22target%22%20variable.) using the scikit-learn python package.
However, when I go to implement my data I get a ValueError where it can not convert string to float, which seems obvious as the tutorial uses two sets of numbers and not strings.
I was wondering if there was a known way to cluster kmers and their counts or simply just strings.

This isn't good data for clustering based on the count value because it needs to be multivariate, e.g. you'd need two or more 'count' columns (if that makes any biological sense). This data is univariate so I'd just use a standard bar chart, personally I would rank the data and consider the underlying distribution.
You could make a tree of the sequence data and the data counts are the labels. If you want to make a tree (hierarchical clustering) then using nucleotide data its becomes a phylogenetic tree, which isn't hierarchical clustering - although clustering can be used to make phylogenetic trees.

K-mers ... if these sites are not homologous a phylogenetic tree is not a good solution: they do look like homologous however. In terms of shifting to numerical variables the technical term is called one-hot-key encoding and scikit-learn has a function for this. This is what Steve is saying. However, understanding what the k-mer represents is kinda important. To perform one-hot-key encoding its:
from sklearn.preprocessing import OneHotEncoder
hot = OneHotEncoder()
.... gets complicated

Alternative approach
import pandas as pd
df = pd.read_csv("kmer.csv")
one_hot_encode = pd.get_dummies(df, columns = ['kmers'])

This code will transform your data to what Steve said .. keeping it as pandas is sort of useful because hierarchical clustering will take a pandas input.
I do have doubts ... the counts column needs a unique identifier in addition to the number, because you'll get identical labels otherwise (my original idea had limitations). Ultimately, I don't precisely know what the kmer value biologically means and thats essential. I would be cautious about dumping that straight into a hierarchical clustering algorithm (it will be in scipy BTW not scikit-learn) .

The following would provide a phylogenetic tree input, but statistically its dodgy without understanding what the k-mer means
import re

kmer = ['TGCCAA', 'TGGACA', 'TGGGAA', 'TGGTAA', 'TGTCAA', 'TGTGAA',  'TTCCAA', 'TTGAAA', 'TTGCAA', 'TTTAAA']

count = [2, 2,  2, 5, 5, 2, 13,  6,  4,  6]

x = [chr(i) for i in range(ord('a'),ord('z')+1)]
n = 0
mydict = {}
for c,k in zip (count, kmer):
    n += 1
    mydict[str(c)+x[n]] = k
    if n > 25:
        n = 0
for key in mydict:
    if mydict[key] == None:
        continue
    print ('>' + key + '\n' + mydict[key] + '\n')

output
>2b
TGCCAA
>2c
TGGACA
>2d
TGGGAA
>5e
TGGTAA
>5f
TGTCAA
>2g
TGTGAA
>13h
TTCCAA
>6i
TTGAAA
>4j
TTGCAA
>6k
TTTAAA

Here's the tree

All I can conclude is count and kmer - with one exception TTTAAA and TTGAAA are not correlated. TTRAAA maybe be a result, possibly.
Here's the matrix.

EDIT
Actually there might be a result here: the k-mer representing '13' is a genetic outlier. Whilst it might be difficult to see in the tree I've presented that could proven more decisively. If a bigger sample size demonstrated comparable phenomena - that is a result for certain. So maybe it's not such a bad analysis after all - it could lead to notably interesting result. Hierarchical clustering BTW couldn't generate that results - its assumption base removes the outlier effect.

I think you need some way to convert the k-mer strings into a numerical representation so that they can be clustered. There's a few ways to do this. For example, you could one-hot (or two-bit) encode each k-mer so that each nucleotide is represented as a set of binary values. Or you could just use a simple integer encoding. For example:
def int_encode(kmer):
    two_bit = {'A': 0, 'C': 1, 'G': 2, 'T': 3}
    encoding = []
    for nucleotide in kmer:
        encoding.append(two_bit[nucleotide])
    return encoding

You could then concatenate the output with the frequency counts of each k-mer. For your set of k-mers, your input for clustering would then look like:
[3, 2, 1, 1, 0, 0, 2]
[3, 2, 2, 0, 1, 0, 2]
[3, 2, 2, 2, 0, 0, 2]
[3, 2, 2, 3, 0, 0, 5]
[3, 2, 3, 1, 0, 0, 5]
[3, 2, 3, 2, 0, 0, 2]
[3, 3, 1, 1, 0, 0, 13]
[3, 3, 2, 0, 0, 0, 6]
[3, 3, 2, 1, 0, 0, 4]
[3, 3, 3, 0, 0, 0, 6]


text:
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

You'll want to use the find command, as the answer from @Ram RS shows. But on a separate note, neither of the above commands are syntactically correct. For a safe way to use these commands & test them, substitute ls for rm, to simply list files rather than delete them.

If you try find with -links and let me know

What operating system are you using? Not all find implementations have the -empty or -delete flags.

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

The potential problem with this approach in Nextflow is there are "extensive" (or more than I think reasonable) numbers of subdirectories in the output and the names get complicated, which precisely is why I avoid it  (I'd rather construct my own Dockerfiles). find will run straight through the directory structure and the depth can be set. Thats why in my opinion it's the favoured solution here.

@M__ yes, absolutely. If the OP's find supports the -empty and -delete, it is the best tool for the job. You could easily make this recursive with shopt -s globstar; for file in  *.fasta */**.fasta; do....

thank you so much for this line of code!! and the comments really help too. I am running this with a mac OS which is why I think the find function doesn't work for me. It definitely seems like the easier option so I might try and find a workaround with an aws instance in the future but for now this works like a charm. Thank you!

@rimo yes, exactly, macOS ships with a different find implementation that might not support all of the flags used like -maxdepth. -empty and -delete.

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

I saw that link too and tried it but it didn't work for some reason

What happened when you tried it?

it doesn't do anything...nothing output in the STDOUT and nothing changed in the directory

Are you running the command from the right directory? If it does not delete anything, but you're sure you have "empty" files, check if the files have blank spaces or some such invisible characters. cat -A them.

@rimo Ram's solution only works on Linux systems, what are you using?

I'm using a mac OS which is probably why that doesn't work for me

It should work on macOS too - it's not because of the OS, it's because of the file type. Empty aliases are not empty files, much like how empty directories are not empty files. Define your tasks more specifically and you'll be able to get to an accurate answer quicker.

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.

@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).

I think the issue I'm facing is that these files are "Alias" files which are the results from a Nextflow script I wrote. Therefore, the file has size (~180 bytes) but is ultimately empty inside. I think the solution I am looking for is more specific than using the "find" command because it's not working for me. I would like to try and use a solution that actively looks in the files and sees whether there is text there (i.e. wc -l < 1) and if not, delete the file.

Or change the scope of what you are doing to only make softlinks to files that are greater than 0 bytes in size, rather than deleting the ones you already have made.

If find doesn't return any results, a pipe to xargs rm -rf could remove everything in your current directory.

@Steve the point was addressed above

if the -delete fails, that is most likely because the OP is using a find version that doesn't have it, so will very likely also not support the -empty. Also, note that this will fail on weird file names (which shouldn't be an issue with fasta files, but you never know). More importantly, your "edit" section does not remove symbolic links, the -links predicate is about the number of hard links. Since most inodes only have a single hard link (the file), your command will simply delete all files in the directory.

@M__ no, it will not remove symbolic links. Or, more accurately, it will remove all files, which includes symbolic links so in that sense it will remove symbolic links. See man find where the -links option is documented and where it explains that it refers to hard links and also just try touch file{1..10} and then run find . -type f -links 1 -print -delete and you will see that all 10 of the files you created will be removed. And yes, the first point applies to all, but the OP commented here hence I left the comment here to explain why it failed for them.

Okay I err ....

I don't understand how an Alias file (or even a symlink for that matter) can be considered empty. Is the alias/link empty when the file pointed to is empty? Or is it empty when the file pointed to does not exist? I'm unable to understand OP's situation.


text:
Trim reads 1kb upstream of sequence
I need a quick way to trim multiple reads in a FASTA file. I need to trim everything that is 1kbp upstream of this sequence AAGAGATGTTCAATCGTTTAAACAAATTCCAAGCTGCTTTAGCTTTGGCCCTTTACTCTCA.
I figure a quick python script might be the way to go but I'm not sure if there's a tool that already does that (i.e. Trimmomatic or something like it) that would be able to do that for me.
Thank you!

What is the source of these reads (i.e. sequencing platform)? Is it necessary to account for some error in the recognition sequence?

Here's one way using biopython. Note that the Seq object has a number of methods that act just like those of a Python string. One of these is the find() method, which returns the index of the first occurrence of the sequence if found. Replace with rfind() if you want the last occurrence if that makes more sense. This solution uses a generator to avoid storing the entire list of trimmed sequences in memory:
import argparse
import contextlib
import gzip
import sys

from Bio import SeqIO

def trim_fasta(fileobj, sequence, num_bases):
    """
    Trim num_bases upstream of the sequence provided
    """
    for record in SeqIO.parse(fileobj, 'fasta'):
        index = record.seq.find(sequence)
        if index != -1 and index > num_bases:
            record = record[index-num_bases:]
        yield record

def open_file(filename, mode='r'):
    if filename.endswith('.gz'):
        return gzip.open(filename, mode)
    else:
        return open(filename, mode)

def get_argument_parser():

    parser = argparse.ArgumentParser(add_help=False)

    group = parser.add_argument_group('input options')
    group.add_argument(
        "fasta",
        type=str,
        metavar="FILE",
        help="The input FASTA file",
    )
    group.add_argument(
        "sequence",
        type=str,
        metavar="STR",
        help="The sequence to find and trim",
    )

    group = parser.add_argument_group('trimming options')
    group.add_argument(
        "-n",
        "--num_bases",
        type=int,
        metavar="INT",
        default=1000,
        help="Trim INT bases upstream of input sequence",
    )

    group = parser.add_argument_group('output options')
    group.add_argument(
        "-o",
        "--output",
        type=str,
        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.output != '-':
            out = stack.enter_context(open_file(args.output, 'wt' if args.force else 'xt'))
        else:
            out = sys.stdout

        input_fasta = stack.enter_context(open_file(args.fasta))
        trimmed_reads = trim_fasta(input_fasta, args.sequence, args.num_bases)

        SeqIO.write(trimmed_reads, out, "fasta-2line")

if __name__ == '__main__':
    main()

usage: trim_reads.py [-n INT] [-o FILE] [-f] [-h] FILE STR

input options:
  FILE                  The input FASTA file
  STR                   The sequence to find and trim

trimming options:
  -n INT, --num_bases INT
                        Trim INT bases upstream of input sequence

output options:
  -o FILE, --output 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


text:
Keeping DNA sequence after changing FASTA header on command line
I have a FASTA header that looks like this:
>7c8250ef-c89f-4d42-9d48-12c8fe245fb2 runid=606f271fc97598006ba5a922136a2c304cef75a5 sampleid=Pool12-1 read=19008 ch=301 start_time=2021-07-03T08:48:18Z barcode=barcode01

And I am able to change it to the desired output here:
>7c8250ef-c89f-4d42-9d48-12c8fe245fb2_001

Using this command:
grep '^>' 001_old.fasta | cut -d ' ' -f 1,8 | sed '/^>/s/$/_001/' > 001_new.fasta

However it completely deletes the DNA sequence for that read and only outputs the names of the headers to the new fasta file:
>7c8250ef-c89f-4d42-9d48-12c8fe245fb2_001
>815440d5-e033-446c-b7dd-5fad4105184a_001
>d36a0fea-101d-4bd9-bee4-adec7fe4d43c_001
>38ca9eb8-bbd2-4d0a-ab06-79d134ad46f1_001
...

I think it has something to do with the cut command I'm using. I tried just doing 1 but that gave the same output so I figured if I included the 8th position of text after a space it would work but again no dice. If someone could help me to include the DNA sequence under the FASTA header in this new file that would be great.

Do yourself a favor and use Bipython for parsing and modifying FASTA files.

The grep, cut, sed pipeline could also be replaced with the following awk command:
awk '/^>/ { print $1 "_001"; next }1' file.fa

The general syntax of an AWK command is condition { action } where condition specifies a condition or regular expression, and action specifies a block of code or set of commands to be run when the condition is met. Both condition and action are optional. Here, we use a regular expression to match lines beginning with a ">" symbol, and print the first whitespace-separated field followed by "_001". We then use next to skip the remaining actions and move to the next line, which avoids any further processing of the current line. The 1 on the end is actually another condition that always evaluates to true. Since there's no action associated with this condition, AWK performs the default action which is to print the entire line. This effectively prints all lines that do not begin with a ">" symbol.

thank you this worked perfectly! and thank you for the explanation too

The problem is in the grep command you're using, which only sends the fasta file headers to the output. Here's my guess at what the commands in your pipeline do:

grep '^>' 001_old.fasta - extract lines starting with '>' from 001_old.fasta

cut -d ' ' -f 1,8 - assume file contains space-delimited records; extract only the first and eighth field from each line

sed '/^>/s/$/_001/' - in lines that begin with a '>', replace the end of the line with '_001'

> 001_new.fasta - write result out to 001_new.fasta

Of all of these, I think only the sed and output parts are important, but it needs to be slightly modified to work as intended:
sed '/^>/s/ .*$/_001/' 001_old.fasta > 001_new.fasta

Assuming I've got the syntax right, this should search for lines beginning with '>', and replace everything from the first space up to the end of the line with '_001'.
I'm more familiar with perl, so here's a perl one-liner that should do the same thing:
perl -pe '/^>/ && s/ .*$/_001/' 001_old.fasta > 001_new.fasta

thank you so much! this makes a lot of sense to me as well.


text:
Convert Symmetrix Matrix into Lower Triangular Matrix
I'm using Scipy's pdist function to calculate a distance matrix which is working great but it gives me a symmetric matrix and I want to convert it into a lower-left triangular matrix.
This is what it prints out for me:
  a   b   c   d
a 0 ### ### ###
b ### 0 ### ###
c ### ### 0 ###
d ### ### ### 0

This is what I want it to be:
###
### ###
### ### ###

I know there's a function in R to do this (upper.tri) but it just replaces everything to the right of the 0 with more 0s. Is there a better/easier way to get it like this? I'm opening the file in excel and removing everything to the right with the 0 manually right now.

From reading the scipy documentation for the pdist function, it seems to suggest it should already only return upper-triangular entries rather than a full matrix. Are you using pdist directly or through some other package?

This is the command I use to calculate the distance: pd.DataFrame(squareform(pdist(kmer_matrix, 'euclidean'))) I convert it into a pandas dataframe so it's easier to output into a file and I use squareform because I read somewhere that it helps space out the final output.

I just tried running it without squareform and it outputs all the distances in a descending list :/

NumPy has a convenient built-in function for accessing the lower-triangular or upper-trianglar entries of an array.
To get the lower-triangular entries, you can use numpy.tril_indices
to get the indices of these entries, and then use these indices to extract the entries.
import numpy as np

ar = np.array(
    [[0, 1, 2],
    [3, 4, 5],
    [6, 7, 8]]
)
idx = np.tril_indices(3, -1)
ar[idx]
# outputs array([3, 6, 7])

This doesn't automatically produce the spacing you want, but formatting the final display is likely context-dependent.

This could work for now thank you!!!

There seems to be a scipy method that retains the lower triangle but with the 0s. See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.tril.html#scipy.sparse.tril

Thanks! This is pretty helpful too


text:
Cell Count Calculation in Split View Using FeaturePlot
In single-cell sequencing analysis, the FeaturePlot function allows us to visualize gene expression levels across different cells, and the split.by parameter enables us to create a split view based on a categorical variable, such as cell clusters. My question is, after employing the split view in FeaturePlot with the command FeaturePlot(seurat_object, features = "gene_of_interest", split.by = "cluster_id"), is it possible to perform cell count calculations for each cluster within the split view?

Yes, it's possible - I've added a cluster/bin biplot a function in my Single Cell Browser application - but graphical region selection is not implemented in any of the Seurat functions.
When I do this myself I prefer to use the command line interface rather than the graphical interface. You can use the table function within R to create count bins for multiple features at the same time:
> table(sc.data$cluster, sc.data$sampleTag)

Annotating Seurat plots (like FeaturePlot) with other graphics is possible, but difficult. I've found it easier to fetch the raw data and create a custom display function.


text:
snRNA-seq Demultiplexing before or after QC
I'm currently analyzing snRNA-seq data and have to demultiplex them using HTODemux in Seurat.
The question is : should I do my QC filters on my droplets before or after the demultiplexing ?
From what I can test, the QC filters actually is going to have an effect on the assignation to each sample, as you can see in the small summary data below
Demultiplex, then QC :
Cutoff for 2-2 : 104 reads
Cutoff for 2-3 : 242 reads
Cutoff for 5-2 : 80 reads
Cutoff for 5-5 : 78 reads

 Doublet Negative      2-2      2-3      5-5      5-2
     347     1754     1182      960     1278      959

QC, then Demultiplex :
Cutoff for 2-2 : 109 reads
Cutoff for 2-3 : 194 reads
Cutoff for 5-2 : 83 reads
Cutoff for 5-5 : 74 reads

 Doublet Negative      2-2      2-3      5-5      5-2
     426     1466     1066     1284     1302      936

As you can see, in the case of sample 2-3 we have a 25% difference in assignation, so I don't think we can completely ignore that.
To be clear, I understand why there's a difference (normalisation, clusterisation differences, quantiles differences), I am looking for best practices here.
Thanks in advance

In the end I took the decision to first do the QC and then Demultiplex, reasonning that the Demultiplex steps (normalisations, clustering etc ...) are probably cleaner this way, leading to a (I hope) better assignation


text:
Naming conventions for kegg nodes with multiple genes
So I'm downloading XML files from Kegg using BioConductor and and I'm running into a problem. If we look at HIF-1 signaling pathway we see the node growth factor (GF) on the left hand side, we can trace it's XML entry to this entry:
    <entry id="17" name="hsa:1950 hsa:3479 hsa:3630" type="gene"
        link="https://www.kegg.jp/dbget-bin/www_bget?hsa:1950+hsa:3479+hsa:3630">
        <graphics name="EGF, HOMG4, URG..." fgcolor="#000000" bgcolor="#BFFFBF"
             type="rectangle" x="87" y="272" width="46" height="17"/>
    </entry>

While there is no GF on this, there are the three info card links that I can use to extract information. However, when I click on that link, there isn't a place on those info cards that says that these three genes are related to that general title of GF.
Normally this isn't a problem, I would just look at the map, but I'm trying to do something at scale where I can't look at every instance of this.
So my question is, how do I find the labels that Kegg uses in its map using the XML file/python package?

Python's Bio.KEGG.Gene has an attribute called name described as "A list of the gene names". Might help? Otherwise "orthology" is its own universe. Last ditch attempt try beautifulsoup and webscraping. In theory webscraping will work because if you can access the links via clicking them beautifulsoup will do this too.

Thanks, I'll check out the name attribute! I'm trying to avoid scrapping at all costs cause this particular use case won't be pretty, but if it's all there is then...

How do you download the files from KEGG with Bioconductor? It might be easier to use Bioconductor to retrieve the gene names

If you used KEGGREST you should be able to retrieve the gene names of each pathway: https://bioconductor.org/packages/release/bioc/vignettes/KEGGREST/inst/doc/KEGGREST-vignette.html

The problem with the name attribute is that it returns multiple names for each card inside that node. GF for example returns insulin, insulin like growth factor 1, and epidermal growth factor. Which one do I use? and how did the original map just label it as GF (growth factor)?


text:
Mass spectrometry: How to obtain expression counts for protein
I was given a data set consisting of 6 samples: 3 control and 3 treated. I dont have information about how data were generated and what the outputs represent.
Each sample corresponds to a folder within which there are the following files:
"DB search psm.csv" "de novo only peptides.csv" "de novo only peptides.xml" "de novo only spectra.mgf" "filter.csv" "peptide features.csv"
"peptide.csv" "peptides_1_1_0.mzid" "peptides.pep.xml" "protein-peptides.csv" "proteins.csv" "proteins.fasta" "PSM ions.csv"
The "peptide.csv" files have the following column names:
Peptide -10lgP Mass Length ppm m/z RT 1/k0 Range Area Ctrl-2 Fraction Precursor Id Source
The "proteins.csv" files have the following column names:
Protein Group Protein ID Accession -10lgP Coverage (%) Coverage (%) Ctrl-1 Area Ctrl-1 #Peptides #Unique #Spec Ctrl-1 PTM Avg. Mass Description
Has anybody ever worked with these type of data?
I would like to do exploit these data and perform a differential analysis but i don't know how to create a SummarizedExperiment from these data.
Any help would be greatly appreciated.
Thank you

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

Could you provide an example dataset? And post the headers and top 5 rows of each file?


text:
Merge/Join one data frame (X) to multiple data frames separately
Let us say that we have some data frames named S1,S2,S3... that I want to join independently by column names to dataframe X. This is how they look like:
S1 <- data.frame(
  values = c (1:5),
  numbers = c(11:15))
S1

S2 <- data.frame(
  values = c (1:5),
  numbers = c(11:15))
S2

S3 <- data.frame(
  values = c (1:5),
  numbers = c(16:20))
S3

#Dataframe X
X <- data.frame(
  values = c (26:30),
  numbers = c(501:505))
X

S1
S2
S3
X

Using rbind, it is really easy to merge/join them according to column names. As below:
S1_X<- rbind(S1,X)

S1_X

S2_X<- rbind(S2,X)

S2_X

S3_X<- rbind(S3,X)

S3_X

Now the problem is that I have 500 plus data frames (S1,S2,S3,S4,S5.....S500) that I want to join to dataframe "X" by any function. For example rbind. Is it possible to use some loop/function for that.
What I have best tried is create a list (nested list of dataframes) that I can later split into individual dataframes:
output = lapply(list(S1=S1,S2=S2,S3=S3),rbind,X)

The output is what I expected. So it works The only problem is, how do I specify S1 to S500?
Thank you in advance!

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

Hi, I did specify that I have 500 data frames named S1,S2,S3,S4.....S500 that I want to join to data frame X which is a dataframe too. If you will run the codes above you will find out that I tried to make it simple. Thank you!

@salman1490 I would recommend deleting this question here and reposting on Stack Overflow. It's not really as bioinformatics question and you'll get a much quicker response there.

You can use For Loop here. Here I have created code using your algorithm, first created a list of merged dataframes and then split the list into separate dataframes.
S1 <- data.frame(values = c (1:5), numbers = c(11:15))
S2 <- data.frame(values = c (1:5), numbers = c(11:15))
S3 <- data.frame(values = c (1:5), numbers = c(16:20))
#Dataframe X
X <- data.frame(values = c (26:30), numbers = c(501:505))

#defining a NULL object
S=NULL
#creating a list of merged dataframes
for (i in 1:3){
S[i] = list(rbind(get(paste("S",i,sep='')),X))
names(S)[i] = paste("S",i,"_X",sep='')}

#extracting dataframes
for(i in 1:3){
assign(paste("S",i,"_X",sep = ""),S[[i]])
}


text:
Avoiding false anticorrelations between individual genes per cell in single-cell RNAseq
When analysing a single-cell RNAseq dataset, one question I like to ask is: At a single-cell level, are these two genes correlated (expressed together) or anticorrelated (only one or the other is expressed)? However, I am concerned that - given the low read count per cell - it may falsely appear that two genes are anticorrelated when in fact it is just unlikely that both RNAs are detected at the same time.
This problem can partially be dealt with using an imputation technique such as SAVER, to estimate the true relative expression values of undetected genes. Additionally, a statistical test (e.g. permutation) can be performed to ask whether the incidence of cells with both transcripts detected is higher or lower than would be expected, given the total number of trancripts.
My concern, however, is that there may be additional technical aspects of the sequencing process that will exaggerate the apparent anticorrelation between genes. For example, in a cell where the total number of RNA transcripts at any one time is low, a single transcript may be read multiple times (due to multiple primers binding, resulting in multiple cDNAs).
In addition to this, are there any other technical or analytical aspects of ssRNAdeq likely to lead to false anti-correlations?


text:
Visualize RNA co-folded secondary structure with Python
I would like to visualize the secondary structure of co-folded RNA strands in Python, preferably in a Jupyter notebook.
What are the recommended tools? Can anyone provide an example?

... . I know both aspects of this question and I'm not aware of an interactive Python package using e.g. iPython magic commands that will do this. The thing you could try is  draw_RNA https://github.com/DasLab/draw_rna
draw_RNA will operate on Jupyter, but I've not tried it.
A very, very good package is VARNA (its Java BTW). The VARNA group make further recommendations from there for complex structures. I would vouch for this package for complex structures.
A good approach is to output in Stockholm 1.0 format. It then just drops and drags into VARNA. Dot-bracket format also works and I would assume ct format works too.
A good reason why in-line graph imaging isn't the best solution for complex structures, i.e. the in-line graphs in a Python debugger or Jupiter is complex RNA structures need to be interactive because you need to move them: they "clump" together and have to manually separate them.
As a final note Forna may be of benefit,


text:
Warning message related to tensorflow package when trying to run deepbgc bioinformatics tool
I have come across some related issues pertaining to this tensorflow warning but none related specifically to deepbgc
I have created a conda environment and have installed onto it the bioinformatics analysis tool deepbgc. Now the program runs but it does output some warning messages related to the tensorflow package. When I run the "deepbgc info" command through terminal, I get the following output messages related to the tensor flow package:
Using TensorFlow backend.
WARNING 03/10 12:55:05   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:74: The name tf.get_default_graph is deprecated. Please use tf.compat.v1.get_default_graph instead.

WARNING 03/10 12:55:05   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:517: The name tf.placeholder is deprecated. Please use tf.compat.v1.placeholder instead.

WARNING 03/10 12:55:05   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:4138: The name tf.random_uniform is deprecated. Please use tf.random.uniform instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:133: The name tf.placeholder_with_default is deprecated. Please use tf.compat.v1.placeholder_with_default instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:3445: calling dropout (from tensorflow.python.ops.nn_ops) with keep_prob is deprecated and will be removed in a future version.
Instructions for updating:
Please use `rate` instead of `keep_prob`. Rate should be set to `rate = 1 - keep_prob`.
WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:174: The name tf.get_default_session is deprecated. Please use tf.compat.v1.get_default_session instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:181: The name tf.ConfigProto is deprecated. Please use tf.compat.v1.ConfigProto instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:186: The name tf.Session is deprecated. Please use tf.compat.v1.Session instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:190: The name tf.global_variables is deprecated. Please use tf.compat.v1.global_variables instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:199: The name tf.is_variable_initialized is deprecated. Please use tf.compat.v1.is_variable_initialized instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:206: The name tf.variables_initializer is deprecated. Please use tf.compat.v1.variables_initializer instead.

Could this be related to the version of tensorflow? Currently, I have version 1.15.4 installed onto the conda environment. Are these error messages significant and could they impact the data running through deepbgc? Or can they mostly be ignored? Any guidance appreciated, thanks!

That should be easy to sort out and its good to get it sorted. I'll reply tomorrow.

@M_ great, I appreciate it. I was thinking maybe this was an issue with my tensor flow version so I tried updating it but the problem with that is it makes other packages like numpy incompatible with the version of deepbgc I'm running on. Could an update to deepbgc do the trick and then subsequently, an update to tensorflow?

The 'blunderbus' way to do this is as follows (not recommended in this instance :-) ),

    conda update --all
    conda update -n base conda
    conda update -n myenv --all

The refined way to do this is

    conda update -n myenv tensorflow # see `Best guess, notes

In theory conda was supposed to sort out the deprecation, because all the versions were supposed to be harmonious. I'm a bit puzzled (see below). What might have happened is tensorflow has been upgraded (or keras) against the version the developer had written it for.
Here myenv is the name you gave that given conda environment.
In your case, I would look at point 2, because the environment is very specific to your project and ensures reproducibility.
Note A purist would specify the exact channel -c that tensorflow is found in, or whatever dependency. This actually can be very important because there can be name conflicts, obviously not for tensorflow because I'm pretty sure Google have the name registered.
The reason I do point 1 is the problem might a dependency of tensorflow and you'd have to figure it out, but again see 'notes'.

Best guess
In this particular example,
conda update -n myenv keras

Might be the ticket, because what I think is going on is its not actually tensorflow but the keras sitting above it.
The alternative (probably better) guess is ....
conda update -n myenv -c bioconda deepbgc

This is probably what I would try first.
Note I'm surprised that conda/mamba has not resolved this issue on installation because thats exactly what is it is supposed to do. In other words its supposed to make a perfect installation of deepbgc precisely as the developer had intended it. There is a risk that updating might not be the solution ... and solution is a new conda/mamba installation (see note 2).
Note 2 Watch for the Python version, you are using 3.7 ... that might be an issue. If the developer had written this for >3.7 that could cause problems. Different Python versions have different methods and if one of the 'new' methods was being used that would be an issue. Ultimately your  package/pipeline is really complex and the more stuff the more chance stuff doesn't fit together. We're currently on 3.10 and 3.9 is really stable.
BTW 3.11 will be a major upgrade of Python3 and will be the 'goto' version.

I looked at deepbgc its built on Python3.7, I don't know why because it was released in 2019. My advice is start again with a new install and follow the instructions exactly as on Github. The order of the build is really important (omitting pip use conda/mamba instead). It doesn't seem to be very conda but thats how the developer did it. In theory the whole pipeline could have been placed on conda. Instructions here https://github.com/Merck/deepbgc
The key issue is that it is built on an old build of Python, with a lot of dependencies which have substancially changed. There is concern about ...
conda create -n deepbgc python=3.7 hmmer prodigal

because the author is not specifying the versions of hmmer or prodigal. These should have been built into the conda upload.

Thanks! I will go ahead and give this a try. Now just want to confirm this, the conda update -n myenv ... command is run outside the environment, correct? In other words, when the environment is deactivated?

Yes outside the given env. In this case I don't think it will make any difference but its you never know.

So when I try running conda update -n my env keras, I am getting the following error: PackageNotInstalledError: Package is not installed in prefix.
prefix: /EFS/tools/miniconda/envs/deepbgc2
package name: keras - should i be providing the full path to the package when i update? I tried to cd into the directory where the package is located and I got the same error message

I've updated the post


text:
KEGG retrieve hierarchy for organism in biopython
I want to obtain a direct graph from KEGG as well as the genes present in a single term. While it is pretty easy to get list of genes/compounds that are part of a specific category, I am kinda stucked to obtain the hierarchical relationships between terms. The docs on biopython.KGML are not really helpful. This is what I have so far
from Bio import SeqIO
import Bio.KEGG.REST as kegg
from Bio.KEGG.KGML import KGML_parser
from Bio.Graphics.KGML_vis import KGMLCanvas
from Bio.Graphics.ColorSpiral import ColorSpiral
from itertools import chain

k =  kegg.kegg_list('pathway', 'pae').read()
todf = []
for x in k.split('\n'):
    if x=='':
        next
    dd = x.split(' - ')
    dd = [v.split('\t') if '\t' in v else v for v in dd]
    dd = list(chain.from_iterable(dd))[:2]
    try:
        ll= kegg.kegg_get(dd[0]).read().split('GENE')[-1]
        ll = [x.lstrip() for x in ll.split('COMPOUND')[0].split('\n')]
        ll= [x.rstrip().split('\t') for x in ll]
        for gn in ll:
            gn.extend(dd)
        todf.append(pd.DataFrame(ll))
    except IndexError:
        pass
todf = pd.concat(todf)
# remove compounds and random strings
todf= todf[todf[0].str.startswith('PA')]
todf[0] = [x.split(';')[0] if ';' in x else x for x in todf[0]]
todf[0] = todf[0].str.split('  ')
todf=todf.explode(0)

and now I have the protein to term relationship but not sure how to get the hierarchy between terms.
Any idea?

Ok realized it was much easier than expected.
On KEGG it is possible to DL a json file, which is intrinsically ordered. So   this parses it into a df
import pandas as pd
import json

todf = []

kegg = json.load(open('meta/KEGG_pae00001.json'))
for main in kegg['children']:
    for broad in main['children']:
        for sub in broad['children']:
            for gn in sub.get('children', [None]):
                if gn:
                     nm = gn['name'].split(';')[0].split(' ')
                     todf.append([nm[0], nm[1], sub['name'], broad['name'], main['name']])
                else:
                     # some terms have no children
                     pass
                     #print(sub)

todf = pd.DataFrame(todf)
todf.columns=['Locus', 'GeneName', 'Cat1', 'Cat2', 'Cat3']
todf.to_csv('meta/kegg_parsed.csv')

Hi, assuming this answer works for you, could you please mark this as the accepted answer? That will help the automatic answer-finding algorithms do the right thing.


text:
How can I assemble my genome from raw files?
I've had my whole genome sequenced (at 30x average coverage) by a lab, and they have provided the raw files to me (BAM, FASTQ, and VCF).
How can I assemble it?
And does assembly provide any further insights that I can't get from the separated fragments?
I'd like to assemble and annotate it, and find insights, but I don't know where to start.

What do you want to achieve by 'assembling' it?  That's usually restricted to random organisms where we don't have a reference sequence. Your VCF/BAM/FASTQ data has already been mapped to an assembled genome so it doesn't really make sense to do an assembly of it. You could certainly annotate it thought, but thats a different question altogether.

It would be useful to know the reference genome.

To add the "do not assemble your genome" answers. There might be specific circumstances when you would actually like to do that. For example, if the 30x you got are HiFi long reads and have a reason to suspect there are some interesting large rearrangements in your genome (rare disease?), it might be informative to actually assemble the genome and compare it to the reference to find these rearrangements.
Of course there is also plenty of possibilities in between (still assuming long reads), for example to map your sequencing data to a pan-genome graph or simply use a mapper that is aware of genomic rearrangements (like ngmlr) and call the structural variants.
But in the end it boils down to two aspects you did not mention in your question: What would you like to learn about your genome?

Given you've tagged this as 'human-genome' there is likely little value in performing genome assembly. De-novo assembly of WGS is generally for non-model organisms without a reference genome available, and would require much higher than 30X coverage (and ideally long reads). By providing the BAM and VCF, they'e already saved you some upstream work unless there are needs to re-process the data from FASTQ

No, assembling won't tell you anything new. When working with organisms whose reference genome is known, all you care about is where your sequences differ from the reference. If you know where your sequences differ, then you can infer the rest, they will be whatever is in the reference.
All you want to do is identify variants, those sites where your genome differs from the reference. That's what the VCF file contains. Bear in mind that the reference isn't some sort of platonic ideal of the perfect human, it's just the collected, pooled DNA of 13 anonymous, apparently healthy individuals. You can expect anywhere from about 3 to 5 million differences in your own genome, the overwhelming majority of which will be completely innocuous.
Since you already have your VCF file, your list of variants, that's all you need to investigate. If you really want to, you can make your assembly by taking the reference genome and then adding the variants in your VCF file, but it won't be of much practical use.


text:
Writing to .cif from BioPandas
My question is twofold. First, I've run into an issue working with .cif files and pandas dataframes. With BioPandas, I can read a .pdb file into a dataframe and then go the opposite direction, writing a .pdb file from a dataframe. With .cif files, I can only write into a dataframe - PandasMmcif does not contain a to_cif function like PandasPdb has a to_pdb function. Are there any alternative methods I could use to get from .cif to pandas dataframe back to .cif?
Second, the reason I am not simply converting all of my .cif files to .pdb is that I am uncertain whether I will lose data in the conversion. Is this fear unfounded? I don't know the specific differences between .pdb and .cif.

You may lose some metadata during mmCIF <-> PDB conversion (what is converted and what is ignored depends on the converter that you use), but it's the same when you load a PDB or mmCIF file into a program, modify it and write it back – some metadata may get lost. The most essential things should be preserved. Unless you have a big mmcif file that can't be converted to a standard PDB due to the limits on the number of chains and atoms in the PDB format (there are workarounds for this used by some converters).

Perfect, thank you!

try Biopyhton https://biopython.org/docs/1.75/api/Bio.PDB.mmcifio.html  class Bio.PDB.mmcifio.MMCIFIO
Bases: Bio.PDB.PDBIO.StructureIO

Write a Structure object or a mmCIF dictionary as a mmCIF file.


text:
select highly variable genes out of dataframe
lets say I have the following dataframe:
Ensembl_ID          A   B   C   D   E   F   G
ENSG00000236601.2   0   0   0   1   0   0   1
ENSG00000237094.12  2   187 442 109 88  144 486
ENSG00000269732.1   0   0   0   0   0   0   0
ENSG00000284733.1   0   0   0   0   0   0   0
ENSG00000233653.3   0   0   0   0   0   0   0
ENSG00000250575.1   1   11  17  5   5   8   17
ENSG00000278757.1   0   0   2   0   1   0   7
ENSG00000230021.9   0   25  30  27  9   8   11
ENSG00000235146.2   0   0   0   0   0   0   0
ENSG00000225972.1   0   7   3016    11  5   14  5
ENSG00000225630.1   1   113 194 76  91  47  94
ENSG00000237973.1   25  1037    9767    1160    886 1321    5220

What I would like to do is only keep the ~top 50 highly variable genes:
I thought it would be like this:  dif_genes <- order(apply(df, 1, var), decreasing=TRUE)[1:50]
but that leaves me with a vector containing variances without genes.

One way is adding a new var column to your original data frame and sorting the updated data frame by the newly introduced column:
df$var <- apply(df, 1, var)
df[order(df$var, decreasing=TRUE),]

the OP wants top 50. I would suggest df[order(apply(df[, -1], 1, var), decreasing=TRUE),][1:50, ] to avoid adding the variance as a sample (in case OP wants to do a heatmap) and removing the names that seem to be in the first column. The key @user16197 to keep the top 50 genes is the comma after the 50; it tells R that this is a 2D data.


text:
How to get the external (cluster) job id and not internal snakemake job id in cluster log files?
Background I am using snakemake and I want to add the external (cluster) jobid to my cluster output and error log files. Below is my shell script.
snakemake --configfiles config/config_resources.yml config/config_parameters.yml --latency-wait 60 --use-conda --printshellcmds --jobname PipeName_{name}.jobid{jobid} --jobs 300 --cluster "bsub -W 480 -n {threads} -M {resources.mem_gb}G -R 'rusage[mem={resources.mem_gb}G]' -o logs/{name}_{jobid}.out -e logs/{name}_{jobid}.err".
But I only get the internal one in de log name from snakemake (e.g. 1 to 5), instead of the one from the cluster.
Question How to I get the external job id in there and not the local one made by snakemake? Or where in my shell script did I go wrong?
Cheers and thanks in advance.

bsub supports some filename placeholders using a percent sign, like %J for the LSF job ID, so you can use that:
snakemake ... --cluster "bsub ... -o logs/{name}_%J.out -e logs/{name}_%J.err"
That way snakemake will format the cluster command string with its wildcards (like name) as it constructs the bsub command, and then bsub will put the job ID in place in those filenames.  Since Snakemake just sees that as one big command to run for the cluster and doesn't know the details other than formatting the text, I think the two shouldn't clash.

Oh my, thank you! That worked!


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

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

I believe bcftools is able to accomplish this.


text:
Liftover variants/vcf called on NCBI reference genome
I have been tackling this problem for a while now, and I hope someone has a better answer.
I need to liftover variants (.vcf) from NCBI hg37 to NCBI hg38. The NCBI remapping tool has been depreciated. For liftover you can use multiple tools, but they all need a chain file, or they only liftover annotation files like gff (like flo). The only chain files available (that I can find) are from UCSC and ENSEMBL, and as far as I know, these reference genomes are different.
I have found some very oblique ways to make chain files here, here, and here. When I finally figured out what was going on, the steps seemed straightforward, but when I try to follow them the amount of additional steps you need to take becomes exponential almost.
At a certain point I also found some perl script from NCBI for remapping or getting the files I need (sorry, can't find the link), which I of course didn't get working because every step I took another dependency failed.
I also tried contacting NCBI, to no avail.
Lifting over is essentially aligning one reference genome to another. Until here I conceptually get it.
I cannot believe this is so hard and I am tired :'-). There must be someone who has already figured this out and knows the easier way to do this. My research should not be about "the way to lift over an NCBI genome". Before I continue this struggle (mainly trying to finish making a chain file as described in the links above), does anyone have anything useful for me? Am I missing something?
Thank you in advance

I don't think ENSEMBL have their own genome version, are you sure?

There seem to be differences between the genome provided by ensembl. It is explained more in the link I provided in my post.

Liftover is hard. While many variants can be lifted over relatively easily, some cannot. First, some variants are not even variants in the other genome (the other genome has the variant allele as the reference). Then you can have things that were assembly errors that have been corrected, and also regions that were simply not covered in the old genome. It just isn't a straightforward problem at all.
A few years ago, in 2017, I tried three separate tools, UCSC liftover, NCBI Remap, and a custom tool built by the company I worked for. I used it to lift over bed regions from hg19 to hg38. I got three separate results. Each of the three tools gave me slightly different output which means that in the best case scenario, one of them was "correct" and two "wrong", but more likely all were "wrong" in different ways.
The only safe way I, personally, have found is to convert variants to protein notation, lift over, and check if the new variant still results in the same protein. That cannot be done for non coding variants, of course.
I found a paper, which I have only glanced through, that compares different liftover tools and suggests a combination of approaches. You might want to have a look: Luu PL, Ong PT, Dinh TP, Clark SJ. Benchmark study comparing liftover tools for genome conversion of epigenome sequencing data. NAR Genom Bioinform. 2020 Aug 6;2(3):lqaa054. doi: 10.1093/nargab/lqaa054.
Even there though, and even after limiting to only ungapped blocks, they couldn't get perfect, 100% liftover:

We showed that 100.0% ungapped blocks were successfully converted (Figure (Figure2B),2B), of which 99.7% were lifted to the correct corresponding region on hg38.

And, while they do seem to have been quite successful, they also point oiut there is no way to get this 100% perfect:

Liftover is a rational solution for conversion between genome assemblies by coordinates when there is a lack of data storage, computing resource, time or human resource as it can reproduce highly accurate results for downstream analysis of epigenome sequencing data. However, it is worth noting that liftover can only assure the accuracy of coordinates and annotation conversion, therefore assumptions of achieving equivalent corresponding sequences genome-wide should not always be made. We clarify the principle underlying this powerful tool and propose a strategy to overcome the limitations of liftover at problematic regions.

WARNING: shameless plug of own software below.
If you have a budget for this, you can use the API made by the company I work for to convert your hg19 variants to p. notation and then query them in hg38 to do the conversion. That should correctly lift your coding variants over, but still won't really help with the non-coding ones. Plus, this is not a free tool, so it is only helpful if you are able to pay to use it.

Thank you @terdon for the elaborate explanation. I indeed understand the problem is hard. I have read many comparison papers and never thought any of the tools accuracy or lifting would be possible for a 100%. I was actually just hoping to be able to liftover something, rather than nothing :'-)

@Dandelion then that paper I linked to should help a lot. There's a nice table with different tools included there which should give you something to work with. If you remember, it would be great if you could come back and post an answer with the solution you ended up choosing. I'd like to know about it too!

Could you use the hg19 --> hg38 chain from here?
As I understand it, the main difference between GrCH37 and hg19 are the format the chromosomes are written (whether there is a "chr" appended in front or not), and that they are mostly interchangeable, although there are of course other small differences. Unless you care about those small differences can you just use the hg19 chain file, and use a homemade script to add/subtract the "chr" from the input/output as needed?

I have used this approach too. However, it is not simply the case of adding "chr" to the names. As you can see from the analysis from GATK (I added the link in the original post) there are differences between sequences as well. Unfortunately, for our specific case, also chromsomes we are particularly interested in.


text:
How to use a .out Kraken output format
I was given some Kraken output files but the file extension is .out. I know how to use .KREPORT files but not the .out file. Is there any way convert it to a KREPORT file or what do you recommend?
Thanks in advance.

The extension of a file does not determine its content. What is the format of your .out files?

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

If I read the file using Notepad++ I can see something like:
C   A00700:50 342 120|120 0:79 821:1 0:36 |:| A:1 0:58 821:4 0:53
U   A00700:50 342 139|106 0:69 226186:4 0:7 226186:5 0:31 |:| A:1 0:25 818:3 0:45
Is it useful?
Thanks.


text:
Identifying highly variable genes in scRNA-seq: Seurat vs M3Drop
Following a single-cell RNA-seq workshop, I created a Seurat object (my_data), normalized the data, and then tried to identify highly variable genes using two different R packages: Seurat and M3Drop.
variable_genes_Seurat <- my_data %>%
  FindVariableFeatures(selection.method = 'vst', nfeatures = 2000) %>%
  VariableFeatures()

variable_genes_M3Drop <- my_data %>%
  GetAssayData('counts') %>% # unnormalized
  NBumiConvertData() %>%
  NBumiFitModel() %>%
  NBumiFeatureSelectionCombinedDrop(ntop = 2000) %>%
  rownames()

I compared the results and found out that the gene lists were pretty different and shared only 588 of 2000 genes:
shared_variable_genes <- intersect(variable_genes_Seurat, variable_genes_M3Drop)
length(shared_variable_genes)
> 588

My questions are:

Why are the results so different?
Which feature list - shared_variable_genes, variable_genes_Seurat or variable_genes_M3Drop - should I use for dimensionality reduction and clustering?

Thank you!


text:
Bulk protein annotation with Prokka
Can prokka be used to annotate a set of protein sequences (so bulk protein annotation), or does it require an entire genome for the annotation?
I've been looking through the documentation but haven't found a clear answer to this. I tested it, and the annotation of my protein sequences worked, but I wanted to check if the program expected a genome and that I can trust the output.
Thanks!

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

Prokka doesn't need whole genomes, I've definitely used it for smaller contigs before so I think you're fine!

Thanks for the reply, that's great!


text:
reliable and popular HLA imputation tools
Could someone recommend me popular and reliable tools for HLA imputation?
I'm starting my master thesis project and I will infer the HLA region of 2000 sample in a case-control study; I searched for benchmarks on PubMed but I really wanted to now if someone experienced in this kind of analysis can help me
ty

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

You at the very least have to specify what kind of data you are working with.

This sounds like a question that is opinion based, without a clear best answer, and an answer that is likely to change in the future. Such questions are difficult to handle  with the stackexchange format. May I suggest you try https://www.biostars.org/ ?


text:
How to add N bases to PyRange object (i.e., shifting PyRanges coordinates)
I would like to add n (e.g., 1k) bases to a PyRange object (i.e., shifting PyRanges coordinates).
How can I do that?
Input:
import pyranges as pr
gr = pr.random(3)

# gr:
# +--------------+-----------+-----------+--------------+
# | Chromosome   |     Start |       End | Strand       |
# | (category)   |   (int32) |   (int32) | (category)   |
# |--------------+-----------+-----------+--------------|
# | chr1         |  10000000 |  20000000 | -            |
# | chr2         |  20000000 |  30000000 | +            |
# | chr22        |  30000000 |  40000000 | -            |
# +--------------+-----------+-----------+--------------+

Desired output:
# magical command[gr + 1000]:
# +--------------+-----------+-----------+--------------+
# | Chromosome   |     Start |       End | Strand       |
# | (category)   |   (int32) |   (int32) | (category)   |
# |--------------+-----------+-----------+--------------|
# | chr1         |  10001000 |  20001000 | -            |
# | chr2         |  20001000 |  30001000 | +            |
# | chr22        |  30001000 |  40001000 | -            |
# +--------------+-----------+-----------+--------------+

You can use the .extend() method. Starting with your data above, gr.extend({"3": -10000, "5": 10000}) would generate the following:
+--------------+-----------+-----------+--------------+
| Chromosome   |     Start |       End | Strand       |
| (category)   |   (int32) |   (int32) | (category)   |
|--------------+-----------+-----------+--------------|
| chr1         |  10010000 |  20010000 | -            |
| chr2         |  19990000 |  29990000 | +            |
| chr22        |  30010000 |  40010000 | -            |
+--------------+-----------+-----------+--------------+

I don't know a way of to apply this "range extension" in a strand specific way (to be able to generate the desired output exactly) but if this is really needed, you should be able to generate 2 gr objects, one for each strand, do the manipulation and bring these together.

OP is not doing anything strand-specific, just adding 1000 to coordinates irrespective of strand. Your answer provides more flexibility than OP needs.


text:
From RNA-seq raw counts to normalized data
I have several RNA-seq datasets. Some of them provide RNA-seq raw counts, some provide FPKM, RPKM and some have transcripts per million (TPM) data. Non of them provide fastq files, all data is processed already. At the end I want all datasets to be normalized to TPM.
I'm using this code in order to normalize raw counts to TPM: (using R)
rpkm <- apply(X = subset(dataset),
                MARGIN = 2,
                FUN = function(x) {
                  10^9 * x / genelength / sum(as.numeric(x))
                })

TPM <- apply(rpkm, 2, function(x) x / sum(as.numeric(x)) * 10^6) %>% as.data.frame()

When RPKM is provided, and no raw counts is available, I use the second line in the same code:
TPM <- apply(rpkm, 2, function(x) x / sum(as.numeric(x)) * 10^6) %>% as.data.frame()

And when FPKM is provided, I use this formula to transform the data to TPM:
TPM = FPKM*X

where X = 1e6/[sum of all FPKM of a sample]
NOTE: genelength is obtained using the biomart package in R, to get the transcript length directly out of ensemble.
However, because I know the steps of TPM normalization in theory, one should firstly normalize to gene length, and then to gene depth. I'm not sure this is what the code does, feeling very skiptical about it even thought this code was given to me by a Phd student.
Can you guys please help me with this question? is my code right or should I alter something?
Thank you!

Regarding the use of TPM, it's a good idea to be aware of the issues associated with it, even if it's your only available option. There are a few different TPM questions that have been asked previously, this answer has a reasonable summary.

However, because I know the steps of TPM normalization in theory, one
should firstly normalize to gene length, and then to gene depth.

No, you should normalize to transcript length.  Which requires you to know the relative abundance of all the transcripts so you can correct for all their lengths.
I do not think it's possible to un-normalize RPKM or TPM without knowing all the lengths involved in doing the calculations.
Let me give an example: here is RSEM output for just one gene.  Note that the software thinks there are three transcripts present for this one gene:
transcript_id  gene_id length  effective_length    expected_count  TPM FPKM    IsoPct
ENST00000373020 ENSG00000000003 3768    3718.54 179.95  13.88   7.94    48.61

ENST00000494424 ENSG00000000003 820 770.54  0.00    0.00    0.00    0.00

ENST00000496771 ENSG00000000003 1025    975.54  11.94   3.51    2.01    12.29

ENST00000612152 ENSG00000000003 3796    3746.54 0.00    0.00    0.00    0.00

ENST00000614008 ENSG00000000003 900 850.54  33.11   11.17   6.39    39.10

ENST00000373031 ENSG00000000005 1205    1155.54 0.00    0.00    0.00    0.00

ENST00000485971 ENSG00000000005 542 492.54  0.00    0.00    0.00    0.00

So, what would you say the effective length should be for the gene level calculations, and would you get that right without having this transcript-level breakdown?

I do have the transcript lengths, as I mentioned above I got them from ensemble.

Do you have transcript level counts?

No.. I assume this is something you get when you have BAM files? Anyway, I don't have that. I was able to get the transcript length which is necessary for the normalization. Do you think the code here provided by Michael Love is better?
https://support.bioconductor.org/p/91218/#91256

No, bam files do not magically give you transcript level counts.  Software like Kallisto and RSEM and Salmon do, and they do all the calculations for you.  Transcript length is not enough.  Genes have multiple transcripts of different sizes, and you can have a mix of isoforms in your samples.  The software I mentioned will magically handle all this for you.  I don't know that it's possible for you to do the magic in reverse.  I'm sure you can wave your hands and do some math and get numbers.  I don't know that those numbers will be accurate.

I see. According to Michael Love, as you said,  there is no way to calculate accurate TPM values out of count matrix. As I mentioned earlier I have no access to BAM files unforunatlely. But in this link I think I got my answer, they discuss this matter, and Michael Love provides here a simple code to normalize TPM out of count matrix. Thank you for your help and explenations. This is the link https://support.bioconductor.org/p/91218/#91256


text:
Reduce Overfit by removing insignificant proteins through PSEA analysis
I'm training a ML model for disease prediction using protein composition as input, but overfit is present. While looking to remove proteins and reduce multi-correlation, as remove composition predictors through PCA have worked but not enough, this problem appeared and managed to use Metascape (1) for it.
I'm using the Protein Set Enrichment Analysis (PSEA, GSEA for genes) for significance of function as is a method to identify classes of genes or proteins that are over-represented in a large set of genes or proteins, and may have an association with disease phenotypes (2). So:

A higher significance coefficient (log10pvalue) between two proteins suggest they have a similar function so one should be removed, or a lower coef shows the need to erase this quasi-independent protein from the analysis?

Is this analysis good enough to reduce overfitting by remove redundant observations?

References:
(1) https://metascape.org/gp/index.html#/main/step1
(2) https://www.sciencedirect.com/science/article/pii/S2376999819300406

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

Overfitting is usually resolved by adding more data and not removing data sets. The technical blurb is 'augmenting' your training data set to revolve overfitting.
It is possible there is a very large multi-gene family involved distorting the classification, thus reducing that could reduce the overfit and cause the algorithm to spread the classification over a much broader group of genes. The general advice is to find an alternative genomic source of proteins and add that, not for this specific question but in reducing overfit per se.
There are technical issues such as:

a parameter split, does the training have this? If not its strongly recommended - depending on the algorithm.
algorithm, for example decision trees are known for overfitting*, whilst naive Bayes would find difficulty overfitting**.

I personally think multiple measures would resolve this. However, general advice is to look at alternative methods and training approaches before looking at slimming down a data set. I would suggest there needs to be good evidence - which would need to be formally presented. You could be 100% correct (probably are), but its not an orthodox approach.
Clearly there is information here that isn't presented. What you've done is a feature selection based on genes and we don't know what that result was nor the gene group that emerged from that analysis - so it is difficult to dissect this precisely. It could be there is a compelling biological reason for concern of overfitting based on the result, but we'll not know that. What I am certain about is the orthodoxy in resolving overfitting is not leveraged here and any data scientist would be aware of this.
*, this is why random forests are considered strong
**, there are a lot of different approaches to ML and how they are classifying is very different, some are more prone to overfitting than others. Naive Bayesian and 'linear regression' are the least sensitive to overfitting.

Please be aware I do ML within the Python framework and its one of the few areas I'm not aware of the ease of cross-over into other platforms. This is simply because its an analysis where Python excels.


text:
The current status of human genome annotation
I was wondering if you can suggest something on the current status of the human genome annotation efforts. I know about ENSEMBL/HAVANA/GENCODE, and I've read the relevant papers; however, I was wondering if RefSeq or UCSC are still actively developed, and are there any competing efforts on T2T assembly, human pan-genome, etc etc?
I know there are quite a few projects that generate both short and long read RNA-seq and other "evidence" from lots of cell lines etc. I guess I was looking for a good and up-to-date perspective or a review type of article.
Thank you in advance!

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

Welcome to the site @AlexanderPredeus. Thanks for your question, there are a number of different questions here and the format we request is focus on a single issue as the basis of the question. You can always ask multiple questions. What we want to avoid is a discussion forum. Please do consider editing your question.


text:
What is the difference between profile and pattern in Prosite?
I found a definition:
"A pattern is thus a qualitative description of a motif in terms of amino acid sequence. The concept of a profile extends this concept, allowing a quantitative description of a motif, by assigning probabilities to the occurrence of a particular amino acid at each position of a motif."
This is example of profile in prosite
Hits for all PROSITE (release 2022_04) motifs on sequence P98073 [UniProtKB/Swiss-Prot]:
Example of profile
PS50024   SEA   SEA domain profile :
54 - 169:   score = 32.979
QSHEARATFKITsGVTYNPNLQDKLSVDFKVLAFDLQQMIDEIFLSSNLKNEYKNSRVLQ
FEN--GSIIVVFDLFFAQWVSD--ENVKEELIQGLEANKssQLVTFHIDLNSVDILDKLT

This is example of pattern in Prosite
PS01209   LDLRA_1   LDL-receptor class A (LDLRA) domain signature :
197 - 221:     [confidence level: (0)]  CIkadlf.CDgevNCpdgsDEDnkm..C

What is the difference in that case? So in profile I have probabilities (in definition I have "profile extends this concept, allowing a quantitative description of a motif, by assigning probabilities to the occurrence of a particular amino acid". I don't see any probabilities in that sequence.
What does it mean "-" in profiles? Maybe any aminoacid or gap?
Why I have sometimes lowercase letter and uppercase letter in profile and pattern? I don't get it. It's a probability?
What is more important profile or pattern?

Hi @MrMorgotheron and welcome. From my point of view its a good question, but a few too many questions, albeit I could be wrong. The "-" definitely means an indel, so you can remove that question to help focus the question more.

I think this was already answered on a different SE site: https://biology.stackexchange.com/questions/46751/difference-between-motifs-domains-patterns-signatures-and-profiles

This publication explains it nicely: https://academic.oup.com/bib/article/3/3/265/239712?login=false


text:
How to count number of motifs in Prosite?
For example I checked motifs in Transferrin P02787  https://prosite.expasy.org/cgi-bin/prosite/ScanView.cgi?scanfile=816050634835.scan.gz
I have
hits by profiles: [2 hits (by 1 profile) on 1 sequence]
and
hits by patterns: [5 hits (by 3 distinct patterns) on 1 sequence]
Can I say that I detect 4 different motifs? One motif is searched by profile searching (PS51408   TRANSFERRIN_LIKE_4   Transferrin-like domain profile two times) and three motifs are searched by pattern searching (PS00205   TRANSFERRIN_LIKE_1   Transferrin-like domain signature 1, PS00206   TRANSFERRIN_LIKE_2   Transferrin-like domain signature 2 two times, PS00207   TRANSFERRIN_LIKE_3   Transferrin-like domain signature 3 two times).
All the motifs found are equal to 7? Is it right? Can I talk about the different motifs found by the profile and pattern, or are they the same in some way?

As per PROSITE database manual:

Unlike patterns, profiles are usually not confined to small regions
with high sequence similarity. Rather they attempt to characterize a
protein family or domain over its entire length. This can lead to
specific problems not arising with PROSITE patterns. With a profile
covering conserved as well as divergent sequence regions, there is a
chance to obtain a significant similarity score even with a partially
incorrect alignment. This possibility is taken into account by our
quality evaluation procedures. In order to be acceptable, a profile
must not only assign high similarity scores to true motif occurrences
and low scores to false matches. In addition, it should correctly
align those residues having analogous functions or structural
properties according to experimental data.
Profiles are supposed to be more sensitive and more robust than
patterns because they provide discriminatory weights not only for the
residues already found at a given position of a motif but also for
those not yet found. The weights for those not yet found are
extrapolated from the observed amino acid compositions using empiric
knowledge about amino acid substitutability.

Given the less biased and somewhat more balanced result provided by a profile I would say that the sequence contains two motifs that correspond to a PS51408 TRANSFERRIN_LIKE_4 Transferrin-like domain profile and provide the bitscore and e-value for the match.
Best,
av


text:
How to run how_are_we_stranded_here?
I'm trying to run the python library called how_are_we_stranded_here. I have installed with pip install how_are_we_stranded_here. I have also installed its dependencies i.e kallisto 0.44.0 and RSeQC. On how_are_we_stranded_here's pypi page, it gives the following command under the usage heading:
 check_strandedness --gtf Yeast.gtf --transcripts Yeast_cdna.fasta \
--reads_1 Sample_A_1.fq.gz --reads_2 Sample_A_2.fq.gz

So I tried running the program in the terminal using my own files, the command was:
    check_strandedness --gtf Glycine_max.Glycine_max_v2.1.55.gtf \
--transcripts Glycine_max.Glycine_max_v2.1.cdna.all.fa \
--reads_1 SRR2079644.fastq

But it keeps giving the error

./how_stranded.sh: line 1: check_strandedness: command not found

I will be very thankful for any help on how to run it.

You might need to include your pip install output too.  The package defines a set of scripts that should get you check_strandedness in your path, and it works for me... so I bet something didn't work right during install.

Would you consider installing this in a conda environment or other virtual environment, to control the background more effectively (with regard to M__'s answer and your comment below). It sounds like there are some issues in the background that might need resolving.

pip should just work, anyway... Two approaches, point 1 is a last resort approach (not recommended), I definitely recommend point 2.
First try,
python3 -m pip install --upgrade how_are_we_stranded_here

Try again does it work? Error logs?, then possibly ...

You need to find out where it is,

   pip show how_are_we_stranded_here
   pip show check_strandedness
   pip show how_stranded.sh

What you really want is the check_strandedness
cd ~
mkdir bin
./pathtoapp/check_stranded # make sure it works
ln -si /pathwaytoapp/how_stranded.sh /home/whatever/bin/
ln -si /pathwaytoapp/check_strandedness ~/bin/check_strandedness

Then write
export PATH="/home/whatever/bin:$PATH"

or modify .profile or .bashrc, then
source ~/.bash_profile

Try again ...

It is better to perform a conda installation

Install either,

Miniconda or;
Anaconda

    conda install -c bioconda how_are_we_stranded_here
    conda create --name stranded
    conda activate stranded
    check_strandedness -h

Thank you so much for taking out the time to give such a thorough reply. It definitely solved the check-strandedness issue but opened a plethora of other issues, so I just gave up.

No problem @pubsurfted, although this one is not 'solved' upvotes always welcome.

upvoted. Thank you!


text:
How to include only the trimmed fastq files into a nextflow channel and exclude any secondary log files created?
I have the following simple nextflow script, that takes fastq paired reads and runs them through the atria adapter trimming.
// Declare syntax version
nextflow.enable.dsl=2

/*
 * pipeline input parameters
 */

params.reads = "$projectDir/data/*_{1,2}.fastq"
params.transcriptome_file = "$projectDir/data/Glycine_max.cdna.all.fasta"
params.gtf_file = "$projectDir/data/Glycine_max.gtf"
params.outdir = "results"

/*
 * define the `TRIM` process that creates trimmed.reads.fastq files
 * given the reads.fastq file.
 */

process TRIM {

    tag "$pair_id"

    input:
    tuple val(pair_id), path(reads)

    output:
    path(pair_id)

    script:
    """
    atria -r ${reads[0]} -R ${reads[1]} -t 1 -l 12 -q 30 -o $pair_id
    """
}

workflow {

    read_pairs_ch = Channel.fromFilePairs( params.reads, checkIfExists:true )

    trim_ch=TRIM(read_pairs_ch)

}

Next what I want to achieve is input the atria.trimmed.fastq files generated by the trim process into a channel so that I could use it as an input for other process. But atria trimming tool also creates two additional log files along with atria.trimmed.fastq. How do i exclude those additional files and only input the trimmed files into a channel for later use.
Best wishes.

Cross posted here: https://stackoverflow.com/questions/74780770/how-to-include-only-the-trimmed-fastq-files-into-a-nextflow-channel-and-exclude

With output: path(pair_id), you're telling Nextflow to expect a file or directory with the value of pair_id in the process working directory when the command exits (successfully). Atria outputs all files to a directory specified using -o PATH1. This directory is created if it doesn't already exist. Note that the output directory defaults to the current working directory if -o is not specified. If you'd like to write the results to a sub-directory, you could try:
process TRIM {

    tag { pair_id }

    publishDir "${params.outdir}/TRIM", mode: 'copy'

    cpus 1

    input:
    tuple val(pair_id), path(reads)

    output:
    tuple val(pair_id), path(trimmed_reads), emit: trimmed_fastqs
    path "${pair_id}/*.atria.log*", emit: logs

    script:
    def (r1, r2) = reads

    trimmed_reads = reads.collect {
        "${pair_id}/${it.baseName}.atria.fastq"
    }

    """
    atria \\
        -r "${r1}" \\
        -R "${r2}" \\
        -o "${pair_id}" \\
        -t ${task.cpus} \\
        -l 12 \\
        -q 30
    """
}

In the above, we used the emit option in the process output definition to assign a name identifier to each declaration. We can use this feature to reference the desired channel in the workflow block. For example:
workflow {

    read_pairs_ch = Channel.fromFilePairs( params.reads, checkIfExists:true )

    trim_out = TRIM(read_pairs_ch)

    trim_out.trimmed_fastqs.view()
}

Results:
$ nextflow run main.nf
N E X T F L O W  ~  version 22.10.0
Launching `main.nf` [distraught_dalembert] DSL2 - revision: aea2d37092
executor >  local (3)
[5f/78908b] process > TRIM (A) [100%] 3 of 3 ✔
[B, [/path/to/work/6a/440a50e9a0725c18f9eede80549cb8/B/B_1.atria.fastq, /path/to/work/6a/440a50e9a0725c18f9eede80549cb8/B/B_2.atria.fastq]]
[A, [/path/to/work/5f/78908b96adee8b8d663d9b4eb4faba/A/A_1.atria.fastq, /path/to/work/5f/78908b96adee8b8d663d9b4eb4faba/A/A_2.atria.fastq]]
[C, [/path/to/work/b0/a3bb3024a4c05fffec555cce5097a7/C/C_1.atria.fastq, /path/to/work/b0/a3bb3024a4c05fffec555cce5097a7/C/C_2.atria.fastq]]

To input the trimmed reads into a downstream process, just ensure that the input tuple matches the input set cardinality declared by the process. If the input declaration is not the same, you may need to apply one or more operators to your channel to make it fit. For your alignment process, you might not need to make any changes. Also, if your alignment process needs to declare one or more reference files, just make sure these are value channels. Most of the time, what you want is one queue channel and one or more value channels when you require multiple input channels, otherwise you might observe unusual behavior (e.g. see this question). The section on multiple input channels in the docs is well worth taking the time to read in my opinion.
A dummy alignment process might look like:
process ALIGN {

    tag { sample }

    debug true

    input:
    tuple val(sample), path(reads)
    path fasta
    path gtf

    script:
    def (r1, r2) = reads

    """
    ls -g "${r1}"
    ls -g "${r2}"
    ls -g "${fasta}"
    ls -g "${gtf}"
    """
}

And the updated workflow:
workflow {

    read_pairs_ch = Channel.fromFilePairs( params.reads, checkIfExists:true )

    transcriptome = file( params.transcriptome_file )
    gtf = file( params.gtf_file )

    TRIM( read_pairs_ch )

    ALIGN( TRIM.out.trimmed_fastqs, transcriptome, gtf )
}

Results:
$ nextflow run main.nf
N E X T F L O W  ~  version 22.10.0
Launching `main.nf` [extravagant_lattes] DSL2 - revision: 29b5315739
executor >  local (6)
[51/e4e261] process > TRIM (A)  [100%] 3 of 3 ✔
[cc/72c53d] process > ALIGN (A) [100%] 3 of 3 ✔
lrwxrwxrwx 1 users 57 Dec 14 00:21 C_1.atria.fastq -> ../../c1/9a85e19c601f2ab2c4af9c6fd203f0/C/C_1.atria.fastq
lrwxrwxrwx 1 users 57 Dec 14 00:21 C_2.atria.fastq -> ../../c1/9a85e19c601f2ab2c4af9c6fd203f0/C/C_2.atria.fastq
lrwxrwxrwx 1 users 40 Dec 14 00:21 Glycine_max.cdna.all.fasta -> ../../../data/Glycine_max.cdna.all.fasta
lrwxrwxrwx 1 users 29 Dec 14 00:21 Glycine_max.gtf -> ../../../data/Glycine_max.gtf

lrwxrwxrwx 1 users 57 Dec 14 00:21 B_1.atria.fastq -> ../../3f/16747fbd2c9bc4c0dc945b24873cca/B/B_1.atria.fastq
lrwxrwxrwx 1 users 57 Dec 14 00:21 B_2.atria.fastq -> ../../3f/16747fbd2c9bc4c0dc945b24873cca/B/B_2.atria.fastq
lrwxrwxrwx 1 users 40 Dec 14 00:21 Glycine_max.cdna.all.fasta -> ../../../data/Glycine_max.cdna.all.fasta
lrwxrwxrwx 1 users 29 Dec 14 00:21 Glycine_max.gtf -> ../../../data/Glycine_max.gtf

lrwxrwxrwx 1 users 57 Dec 14 00:21 A_1.atria.fastq -> ../../51/e4e261d019c65377d63826963b19b1/A/A_1.atria.fastq
lrwxrwxrwx 1 users 57 Dec 14 00:21 A_2.atria.fastq -> ../../51/e4e261d019c65377d63826963b19b1/A/A_2.atria.fastq
lrwxrwxrwx 1 users 40 Dec 14 00:21 Glycine_max.cdna.all.fasta -> ../../../data/Glycine_max.cdna.all.fasta
lrwxrwxrwx 1 users 29 Dec 14 00:21 Glycine_max.gtf -> ../../../data/Glycine_max.gtf

I appreciate your reply so much and I am very thankful for it. I wanted to confirm one more thing as I just started using nextflow. How do I input these trimmed_fastqs into a new process, say, for alignment? Do I use trim_out.trimmed_fastqs channel as input for the alignment process?

@pubsurfted Correct. I've tried to give an example above of what you might need. Please let me know if I've missed anything, thanks!