text:
Conversion FASTA to BED fails
I want to use conform-gt to get rid of some REF prefix differences see here. Conform-gt needs the reference FASTA file split by chromosome and in .vcf.gz format. How can I do this?
I already tried faidx with
faidx --transform bed reference/GRCh38_full_analysis_set_plus_decoy_hla.fa > reference/GRCh38_full_analysis_set_plus_decoy_hla.bed

but this produced the following error:
Traceback (most recent call last):
  File "/home/oem/.local/bin/faidx", line 8, in <module>
    sys.exit(main())
  File "/home/oem/.local/lib/python3.9/site-packages/pyfaidx/cli.py", line 202, in main
    write_sequence(args)
  File "/home/oem/.local/lib/python3.9/site-packages/pyfaidx/cli.py", line 51, in write_sequence
    outfile.write(transform_sequence(args, fasta, name, start, end))
  File "/home/oem/.local/lib/python3.9/site-packages/pyfaidx/cli.py", line 121, in transform_sequence
    line_len = fasta.faidx.index[name].lenc
KeyError: 'HLA-A*01'

Can you please [edit] your question and explain what you mean a bit more? Your title talks about converting fasta (sequences) to bed (regions) which is odd, but could be done if you mean "read the fasta and tell me what regions of a reference genome they map to". But then your question talks about converting fasta to "vcf.gz" format and this doesn't make sense unless you want some sort of gVCF file with an entry for every position of a genome.

I reconsidered my question. I don't think that I really fully understood what I want to achieve with this format-conversion nor if I really understood what conform-gt really does. Back to the drawing board I guess. I will close this question and hope that the original question in my OP can solve my problem.


text:
Filter with bcftools
I want to filter a SNP, specifically CHROM1:POS:630128. I tried to use bcftools for this.
I came up with the following variant of a bcftools call:
bcftools filter -e "CHROM=1&POS=63018" /home/oem/Documents/GitHub/Genomics_prac_guide/trio/files/genome_file.vcf.gz

which leads to:
Error: cannot use arithmetic operators to compare strings and numbers

and the call with unquoted parameters leads to
$ bcftools filter -e CHROM=1&POS=63018 /home/oem/Documents/GitHub/Genomics_prac_guide/trio/files/genome_file.vcf.gz

sh: 1: /home/oem/Documents/GitHub/Genomics_prac_guide/trio/file/genome_file.vcf.gz: Permission denied
Failed to read from standard input: unknown file type

Q1: Is bfctools filter the correct tool here to use?
Q2: Regarding this comment the first error might indicate that I have to annotate the CHROM & POS Tags in the header file. Is this correct?

An example of the file:
##fileformat=VCFv4.2
##fileDate=20221104
##source=PLINKv1.90
##contig=<ID=1,length=248928329>
##contig=<ID=2,length=242157857>
##contig=<ID=3,length=198157342>
##contig=<ID=4,length=190106194>
##contig=<ID=5,length=181288140>
##contig=<ID=6,length=170595721>
##contig=<ID=7,length=159331484>
##contig=<ID=8,length=145075470>
##contig=<ID=9,length=138211490>
##contig=<ID=10,length=133659511>
##contig=<ID=11,length=135075227>
##contig=<ID=12,length=133261768>
##contig=<ID=13,length=114341522>
##contig=<ID=14,length=106873816>
##contig=<ID=15,length=101922277>
##contig=<ID=16,length=90172890>
##contig=<ID=17,length=83203771>
##contig=<ID=18,length=80252738>
##contig=<ID=19,length=58586567>
##contig=<ID=20,length=64323573>
##contig=<ID=21,length=46679699>
##contig=<ID=22,length=50776369>
##INFO=<ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on real reference genome">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  FAM001_ID001
1   69869   rs548049170 T   .   .   .   PR  GT  0/0
1   74792   rs13328684  N   .   .   .   PR  GT  ./.
1   630128  rs9283150   A   .   .   .   PR  GT  0/0
1   791532  i713426 N   .   .   .   PR  GT  ./.
1   792461  rs116587930 G   .   .   .   PR  GT  0/0
1   817341  rs3131972   G   .   .   .   PR  GT  0/0
1   818725  rs12184325  C   .   .   .   PR  GT  0/0
1   820888  rs12567639  A   .   .   .   PR  GT  0/0
1   823656  rs114525117 G   .   .   .   PR  GT  0/0
1   841166  rs12124819  G   A   .   .   PR  GT  0/1
1   858952  rs12127425  G   .   .   .   PR  GT  0/0
1   866156  rs79373928  T   .   .   .   PR  GT  0/0
1   880041  rs72888853  N   .   .   .   PR  GT  ./.

Can you give us an example of your file?

Done! See edited post.

First of all, bcftools might be a bit of overkill if you want one specific variant. All you need is GNU grep (the default on Linux):
zgrep -P '^(#|1\t630128\t)' file.vcf.gz

Or, with any grep:
zgrep -E '^(#|1'$'\t''630128'$'\t'')' file.vcf.gz

That will print all lines that start with a # (the headers) or that start with 1, followed by a tab, then 630128 and another tab.
Alternatively, you can use awk:
zcat file.vcf.gz |
    awk -F'\t' '/^#/ || ( $1==1 && $2==630128)'

All that said, yes, you can also do it with bcftools and in most cases that is the right tool for the job so you might as well use it here too to get used to it. However, you need the right command and that's view, not filter:
bcftools view -r 1:630128 /home/oem/Documents/GitHub/Genomics_prac_guide/trio/files/genome_file.vcf.gz

Note that the file needs to be indexed, so if you have not done so already, run
tabix file.vcf.gz

Also, the reason your unquoted command failed as it did is because the & has special meaning to the shell, it means "run this command in the background". Therefore, when you ran bcftools filter -e CHROM=1&POS=63018 file.vcf.gz that actually ran two separate commands:
bcftools filter -e CHROM=1&

And then:
POS=63018 file.vcf.gz

Which means "set the variable POS to the value 63018 and execute the command file.vcf.gz. Since the file isn't executable, you got:
sh: 1: file.vcf.gz: Permission denied

First of all thanks for the quite detailed answer! This gives me a good introduction into bcftools environment. So when I use view to filter SNP's the reason for  -filter is then to apply some standard filter to the file?
Second question: If I want to remove this single variant from the file is the correct way to do this bcftools view -r 1:630128 -e file.vcf.gz ?

I think filter is more to get things like "show me variants where the AB is X and the GT is Y" and things like this, while in this case you are just focusing on a specific region, albeit a region of length 1, so view is the command for this. Reversing the match is a completely different question. I don't know if bcftools can do that, but zcat file.vcf.gz | grep -vP '^1\t630128\t' will get rid of all variants not matching that specific position.

source: Stack Exchange

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

sort: votes


text:
Renaming samples in vcf file
I am currently asking myself how to rename samples in a vcf file. This is necessary, because the naming of several samples included in thr vcf file was unfortunately not be chosen by me.
I am loading the vcf file via readVcf() from the VariantAnnotation package in R.
To check the sample names, I can use samples(header(vcf.file)), which will give me the names as characters.
It looks like this:

and so on..
(I am actually very sorry that I cannot give a reproducible example here :( )
I also have a vector with all the new names. For example
new.names <- c("1","2","3","4")

I tried to exchange the "old names" with the new.names vector:
samples(header(vcf.file)) <- new.names

which leads to the error
Error in samples(header(data.vcf)) <- names.sample :
  could not find function "samples<-"

I also tried to use the rutilstimflutre package, which has a already designed function for this, but unfortunately the package is outdated and not available for the newest R version.
Is there a way to do this? Perhaps a better way to assign the new names to the old ones?
Thanks a lot.

Easiest to use bcftools, if samples.txt has the new sample names:
bcftools reheader -s samples.txt -o new.bcf old.bcf

Thank you, since i only have a brief knowledge of R, is there an R solution?

Maybe - I'm not sure. This approach doesn't require any knowledge of any programming - just doing this should work perfectly. Generally I prefer not to use R for working with .vcfs because it's slow to load them into memory. If you are working with .vcfs I would highly recommend getting familiar with bcftools as it is very fast and memory efficient.

This proved to be so simple that it took a while to discover: Just use colnames.
x = readVcf("in.vcf")
colnames(x) <- c("new1", "new2")
# At this point, colnames(x) and samples(header(x)) will show different names

# ... but if you write it out, it should use the colnames for both:
writeVcf(x, "out.vcf")
y = readVcf("out.vcf")

# samples(header(y)) and colnames(y) should both be "new1","new2" here

Thank you for posting a solution in R.

Here is an example in python how to rename samples using pyvcf submodule I wrote:
>>> from fuc import pyvcf
>>> data = {
...     'CHROM': ['chr1', 'chr2'],
...     'POS': [100, 101],
...     'ID': ['.', '.'],
...     'REF': ['G', 'T'],
...     'ALT': ['A', 'C'],
...     'QUAL': ['.', '.'],
...     'FILTER': ['.', '.'],
...     'INFO': ['.', '.'],
...     'FORMAT': ['GT', 'GT'],
...     'Steven': ['0/1', '1/1'],
...     'Sara': ['0/0', '0/1'],
... }
>>> vf = pyvcf.VcfFrame.from_dict([], data)
>>> # vf = pyvcf.VcfFrame.from_file('your_file.vcf')
>>> vf.df
  CHROM  POS ID REF ALT QUAL FILTER INFO FORMAT Steven Sara
0  chr1  100  .   G   A    .      .    .     GT    0/1  0/0
1  chr2  101  .   T   C    .      .    .     GT    1/1  0/1
>>> vf.df = vf.df.rename(columns={'Steven': 'John', 'Sara': 'Rachel'})
>>> vf.df
  CHROM  POS ID REF ALT QUAL FILTER INFO FORMAT John Rachel
0  chr1  100  .   G   A    .      .    .     GT  0/1    0/0
1  chr2  101  .   T   C    .      .    .     GT  1/1    0/1
>>> # vf.to_file('renamed.vcf')

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


text:
OpenAPI-Generator Python Client Instructions Examples
Trying to learn a little bit of Python and Biopython I ended up figuring out what a Rest API is when the RCSB PDB repository released its new Rest API.
I tried to retrieve info from the PDB repository using python3 requests,
then I managed to use OpenAPI-generator to generate what they call a 'Python3 Client'.
Now I am at loss I have this ApiClient class:
print dir():

 ['NATIVE_TYPES_MAPPING', 'PRIMITIVE_TYPES', '_ApiClient__call_api',
  '_ApiClient__deserialize', '_ApiClient__deserialize_date',
  '_ApiClient__deserialize_datetime', '_ApiClient__deserialize_file',
  '_ApiClient__deserialize_model', '_ApiClient__deserialize_object',
  '_ApiClient__deserialize_primitive', '__class__', '__delattr__',
  '__dict__', '__dir__', '__doc__', '__enter__', '__eq__', '__exit__',
  '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__',
  '__init__', '__init_subclass__', '__le__', '__lt__', '__module__',
  '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__',
  '__setattr__', '__sizeof__', '__str__', '__subclasshook__',
  '__weakref__', '_pool', 'call_api', 'client_side_validation',
  'close', 'configuration', 'cookie', 'default_headers',
  'deserialize', 'files_parameters', 'parameters_to_tuples',
  'pool', 'pool_threads', 'request', 'rest_client',
  'sanitize_for_serialization', 'select_header_accept',
  'select_header_content_type', 'set_default_header',
  'update_params_for_auth', 'user_agent']

and don't know how to use it. Googling for examples or searching for info on the OpenAPI-generator website doesn't seem to work for me, I find plenty of tutorials and examples on how to generate the client but none about how to use it.
The generated client class should help to better interact with the server as I understood (I could be wrong here). I am asking for helps and hints here because the OpenAPI webpage under Community points to stackoverflow under [openapi-generator] tag. Any tutorial, easy documentation that could help me understand the output of the generator?

OK fine, found useful information on READme.md generated by OpenAPI-generator itself
and tried pdoc to get html document of my OpenAPI-generator generated client Module
Moreover on stackoverflow How to generate API Documentation using openapi-generator? you can find how OpenApi Generator can do generation of html documentation from your yaml or json definition of your API.
For openAPI json api specification I found RapiPDF too : Web Component to generate PDF from OpenAPI Spec (json yml file).


text:
DESEQ2:Error in rownames, what is the problem?
I am using DESEQ2 for RNAseq analysis but i do not understand why I get this error.
> dds <- DESeqDataSetFromMatrix(countData = data,
+                               colData = labeldata,
+                               design =~ subtype, tidy = TRUE)

**Error** in `.rowNamesDF<-`(x, value = value) :
  duplicate 'row.names' are not allowed

I really do not know what is the issue. Could you please help?
Thanks

Hi @athanasia and thansk for your post. Members here do experience of DESEQ2 for transformation (I don't) and whilst the information could be enough, it is etiquette to provide as much detail as you can. For example, your data, size dimension, runs etc. ...

Please edit your question to include the output of dput(head(data)). thanks

The hint is here:
duplicate 'row.names' are not allowed
rownames(data) should return a vector of only unique values, since it doesn't make sense to have the same gene/transcript in two different rows of your count matrix. How did you generate it?
What does which(duplicated(rownames(data))) return? Maybe that'll give you a hint on what went wrong. In general, if you use salmon or kallisto for counting, you should use the excellent tximport package to get the counts into R, then create a DESeq dataset with DESeqDataSetFromTximport().
EDITS:
The replies to my answer prompted me to re-read your post and I saw that you had set tidy=TRUE, which means the duplicates are actually in your first column of your data.

I thought along the same lines, but data cannot have duplicate row.names - no object can have duplicate row.names, so this is some generated object that is being assigned duplicate row names owing to something wrong with data or labelData that has nothing to do with the row names of either object.

Of course, they can. Here's a oneliner for you: m <- matrix(0, nrow=3, ncol=3, dimnames=list(c('a', 'a', 'a'), c('b', 'b', 'b'))). Anyway, I tried and DESeq2 handles this case semi-gracefully by appending numbers to the dupes. However, @Athanasia had set tidy=TRUE, which prompts DESeq2 to look for row names in the first column of the input, where duplicates trigger an error. I can post an MRE for this exact behavior if you'd like.

Thank you, TIL.


text:
PCA plot in R coloured by sample type
I'm a biologist, not a programmer so please be gentle.
So I have a dataset that looks like
Genes  Patient1   Patient2   Patient3
A          324      433         343
B          431       342        124
Z          232       234        267

then I have the sample sheet where it contains sample info like:

Patient1 - Healthy
Patient2 - Disease
Patient3 - Healthy
(110 samples)

I am using:
library(ggfortify)
df <- dataset
pca_res <- prcomp(df, scale. = TRUE)

autoplot(pca_res)

Then I want to do
autoplot(pca_res, data = ?, colour = '?')

I wish to use the info from the sample sheet to colour my PCA based on the state (healthy/disease) using the autoplot function. Is there a way to do this?

I think you need to provide more information for a minimal reproducible example.

I've found this the way to work with autoplot (ggfortify's prcomp autoplot) - Use a df that has the field you need for color as the first column (to make things easy), then do:
autoplot(prcomp(df[,-1], data=df, color='column1_name')

If the color-by field is not the first column, do:
autoplot(prcomp(df[,!(colnames(df) == "color_by_column_name")], data = df, color = "color_by_column_name")

EDIT:
That being said, you're going to need to process your data so patients are along rows instead of column. Transpose your dataset so rows are columns and columns become rows; change patient ID to rownames and Gene name to colnames, then merge the patient category dataset with the expression dataset. Once you merge, your dataset should look like:
         Category    A      B      Z
Patient1 Healthy     324    431    232
Patient2 Disease     433    342    124
Patient3 Healthy     343    124    267

And then,
autoplot(prcomp(df[,!(colnames(df) == "Category")]), data = df, color = "Category")

Hi Ram,Thanks a lot for taking the time to answer me. I did exactly what you told me above and there is an error
Error in colMeans(x, na.rm = TRUE) : 'x' must be numeric
Any ideas why is that ? I am not sure at all

Ensure that the columns you're passing to prcomp are all numeric. That's probably what's going wrong here.

sorted, thanks a lot :)

Can't say I've ever used autoplot, but this is fairly easy to achieve using "base" ggplot.
One thing autoplot does is put the axis on the same scale reflecting the percentage variation captured by each PC - I think this is a good thing, but most PCA plots you'll see don't do this.
As you haven't given us any info on your metadata, I'll just make a few assumptions - lets make it a data.frame called "colData" with a "Status" column.
prcomp stores the PC scores in "pca_res$x"
# generate some metadata
colData <- data.frame(row.names=rownames(pca_res$x),status=as.factor(sample(c("diseased","healthy"),nrow(pca_res$x),replace = T)))

# get % variation for each PC
pca_res$percentVar <- pca_res$sdev^2/sum(pca_res$sdev^2)

# change scale to reflect % variation
d <- t(data.frame(t(pca_res$x)*pca_res$percentVar))

# merge metadata and modified PC scores by row.names (by=0)
d <- merge(d,colData,by=0)

# make the base plot (PC1 vs PC2)
g <- ggplot(d,aes(x=PC1,y=PC2,colour=status))
# add the points
g <- g + geom_point(size=2)
# add a "wow" factor to the points
geom_point(size=4,alpha=0.5)
# display the plot
g

There may be a way to do this with autoplot - e.g.
autoplot(pca_res,colour=as.numeric(colData$status))

gets you part of the way there, but I think using the standard ggplot syntax is easy enough anyway.


text:
How can atomic RMS differences / atomic RMS distributions help us understand the structure of a protein?
For example, I have a paper that says: "the atomic RMS distribution about the mean coordinate positions is 0.41 +/- 0.04 A for the backbone atoms and 0.82 +/- 0.04 A for all atoms"
And I have no idea what any of this means, I found several papers with this kind of information but none of them explain why this is relevant in order to understand the structure of a protein.

Hi @Ig009 there is nothing wrong with the basis of your question. The resolution in question is excellent if correctly reported. Could you elaborate on the basis of your question/investigation please?

The "RMS distribution about the mean coordinate positions" could be restated as "certainty of the structure". In other words, the main backbone that defines the large-scale shape in the structure is expected to be accurate to within about 0.41 angstrom, while the atoms more on the periphery (or otherwise not forming the backbone of the structure) are less accurately solved. How important this is will depend on what you want to do with the structure. If you want to get an idea for what it looks like then the accuracy is less important, if you want to do docking then the more accurate the atoms around any sort of pocket/interaction surface are the better.


text:
How to remove low abundance and less prevalent data from my dataset?
I have a relative abundance table in .tsv format where samples are in column and rows contain the features (pathways). Something like this reproducible example:
# Pathway   ERR2601 ERR26014    ERR260145_Abu   ERR260  ERR2601 ERR2601 ERR2607 ERR2600 ERR26
1CMET2-PWY: N10-formyl-tetrahydrofolate biosynthesis    0.00826631  0.00946472  0.00962571  0.00411153  0.00811879  0.00844421  0.00915602  0.00789428  0.0091798
3-HYDROXYPHENYLACETATE-DEGRADATION-PWY: 4-hydroxyphenylacetate degradation  0   0   0   0.000695608 0   0   0   0   0
7ALPHADEHYDROX-PWY: cholate degradation (bacteria, anaerobic)   0   0   0   0   0   0   0   0   0
AEROBACTINSYN-PWY: aerobactin biosynthesis  0   0   0   0   0   0   0   0   0
ALLANTOINDEG-PWY: superpathway of allantoin degradation in yeast    0.0000331927    0   0.0000335389    0   0   0   0   0.0000189167    0.0000908934
ANAEROFRUCAT-PWY: homolactic fermentation   0.00682776  0.00681683  0.00675732  0.00495716  0.00691969  0.0053589   0.00514634  0.00779249  0.00426403
ANAGLYCOLYSIS-PWY: glycolysis III (from glucose)    0.00866693  0.0100437   0.00829714  0.00524566  0.00960262  0.0075093   0.00716347  0.00916205  0.00866289
ARG+POLYAMINE-SYN: superpathway of arginine and polyamine biosynthesis  0.00106395  0.000544088 0.000774645 0.00370237  0.00102368  0.00133848  0.00148691  0.00114062  0.00188739
ARGDEG-PWY: superpathway of L-arginine, putrescine, and 4-aminobutanoate degradation    0   0   0   0.00164625  0.0000278133    0   0   0   0
ARGININE-SYN4-PWY: L-ornithine de novo biosynthesis 0.00166327  0.00383993  0.000105528 0.00230283  0.00404469  0.00115975  0.000460711 0.00101816  0.00367102
ARGORNPROST-PWY: arginine, ornithine and proline interconversion    0.000132331 0.000244069 0.0000796126    0.000479006 0.000237704 0.000290783 0   0.00066458  0.0000581229
ARGSYN-PWY: L-arginine biosynthesis I (via L-ornithine) 0.0124159   0.00955758  0.0118644   0.00533077  0.00774959  0.0098352   0.0127232   0.00996349  0.0108839

Now, I want to keep features (i.e. pathways) which are present with abundance >0.0001 and present in at least 10% of samples. What bash command can I use in order to do that?

Your "reproducible example" is not reproducible. It should be a tab or comma separated plain text content pasted directly in the post and formatted as code (or uploaded on a GitHub gist). A screenshot of an excel spreadsheet is quite far from reproducibility.

Thanks Ram RS for your response. I shared a google drive link. I didn't know it will be shown as an image. However, I have edited the post and entered the table as code. But, still I don't know why it is not showing as a  .csv table. Please tell me whether you can help me. Thanks for your time.

Your data is a little difficult to understand. Can you describe it better please? Please explain what each row and each column represents. I see duplicate column names (ERR2601) as well as inconsistent prefixes (such as the _Abu prefix on one which seems to indicate some sort of difference from the rest of the columns). If you're looking to pick all rows where at least 10% of columns 2-10 are >0.0001, you should be able to do that in R by counting number of columns that are >0.0001.

Column names are just subject IDs and rownames are the pathway names. So, please don't think about them. Each column denotes the relative abundance of respective pathways. If you add the values of a column, it will give you 1.00. Now, here, as example, I have 9 columns, right? Now, 10% of 9 is 1 (apprx.). So, I want to get an output file containing rows with at least one value > 0.0001.  That means, rows with 1 or more values > 0.0001 will be kept and rows which does not have at least one value >0.0001 will be removed. As, there is not much complex calculation, I would prefer to do in bash.

Bash is the wrong tool for this.

I would look at Python pandas or an R dataframe. The best bash solution would require awk but you'd need to know how to program with it.

Here is a bash/awk solution:
awk -F'\t' 'NR==1{print}; NR>1{N=0;for(i=2; i<(NF+1); i++){if($i >= 0.0001) N++}; if(N>0.1*NF) print}' input.tsv

We make sure the input field separator is tab, then we loop on the number of fields (NF) and add to the tally the N of fields above threshold. At the end we simply print the entire line if N is high enough.
EDITED: for the header, we print the first line (NR==1) as is, then apply the expression when Number of Records > 1.

Thanks a lot @Alexlok. It's working. But, it is also removing the header (i.e. the first line of the table i provided). Is there a way to retain it in the output?

Indeed, my tail command was excluding the header, I edited to apply a different awk rule for the header or the rest of the text.

Thanks again @Alexlok. It's working really good!!!

Assuming you have this stored in a matrix called m, where the pathway labels are row names:
IDX = apply(m, 1, function(x) sum(x > 0.0001)/length(x) >= 0.1)
m_subset = m[which(IDX),]

OP isn't using R. OP also has duplicate colnames.

OP needs to use the right tool for the job, it doesn't matter if they prefer bash.

I agree with this. You'd need to be an @terdon standard coder to be doing this routinely and he'd likely prefer to use awk, in this example. Few coders are that standard. The OP could use awk if they wanted a bash solution.


text:
How can I generate ASV file from nanopore sequencing data?
I am converting the fast5 file to fastq by using guppy basecaller after that by using kraken2 classified the sequence.
Now I am trying to generate ASV file.
Is this possible to generate ASV file from kraken output?
Or is this any other way to generate ASV file from nanopore sequencing?
I am new to nanopore sequencing data. I need help or suggestions to generate ASV file.
Thank you

What is an ASV file? You should add information about that file format to your question. The first information I get talks about MATLAB, which doesn't seem correct.

ASV: amplicon sequence variant

Dear @twinkle please do check the answer below and do consider upvoting/accepting. Its a good answer and it does help with site stats.

I am not sure that kraken/kraken2 does anything similar to this.
Usually to get OTUs/ESVs I think that you want to use a tool that clusters your sequences into groupings that represent the same sequence, allowing for technical error.
Examples might include dada2 or dadasnake, or qiime. This is after the technical steps for adapter/trimming/denoising steps.
kraken2 is more for classifying into known taxonomic categories based on a reference sequence, which is a related but distinct task.
I don't think that I understand what the ESV file format is supposed to look like. This resource seems to have some info on processing data into ESVs, but it is not for nanopore. Possibly it is still useful for the final step of creating that file format. A linked repository has a ton of one-off scripts that may be helpful but are not very well documented. For example, here is a script that they use to generate TSV files which may be the "ASV"/"ESV" format:
#!/bin/bash
source activate qiime2-2018.8

for item in `ls */*table.qza`
    do
    name=`echo $item | cut -d\/ -f1`
    qiime tools export --input-path $item --output-path EXPORTS-GTDB-tax/$name
    mv EXPORTS-GTDB-tax/$name/feature-table.biom EXPORTS-GTDB-tax/$name.biom
    rmdir EXPORTS-GTDB-tax/$name
done

qiime tools export --input-path 11_ANT-28-5-SPA_seqs_classified_against_GTDB_r86.1_807F-1050R_primer_sliced/classification.qza --output-path EXPORTS-GTDB-tax/

sed -i '1c#OTUID    taxonomy    confidence' EXPORTS-GTDB-tax/taxonomy.tsv

#Tutorial here: https://forum.qiime2.org/t/exporting-and-modifying-biom-tables-e-g-adding-taxonomy-annotations/3630
for item in `ls EXPORTS-GTDB-tax/*biom`
    do
    filestem=`basename $item .biom`
    biom add-metadata -i $item -o EXPORTS-GTDB-tax/$filestem.with-tax.biom --observation-metadata-fp EXPORTS-GTDB-tax/taxonomy.tsv --sc-separated taxonomy
done

for item in `ls EXPORTS-GTDB-tax/*.with-tax.biom`
    do
    filestem=`basename $item .biom`
    biom convert -i $item -o EXPORTS-GTDB-tax/$filestem.tsv --to-tsv --header-key taxonomy
done


text:
How can I extract sequences from FASTA file using OTUs ID?
I have one FASTA file with sequences and headers, and one OTUs table.
What I want to do is replace the current headers in the FASTA file with just the OTU IDs from the .txt file. The files look like something like this:
FASTA file:
>9d4544fa-322a-4fac-bcf5-f8ca410c7e8a runid=b48dcb2029ef589b8cdb55fe529d93b3baadbc7c sampleid=28012021COVID3 read=57150 ch=310 start_time=2021-01-28T18:23:35Z
TCAGTTACGTATTGCTGGTGCTGAAGAAAGTTGTCGGTGTCTTTGTGTTAACCTTTCTGTTGGTGCTGATATTGCAGTCGCCCCCACCTCCTTGGATCCGCCAGGTTAAACACCCAAGCAGACGCCGAAGATAGAGCGACAGGCAAGTAGGTTAACACAAAGACACCGACAACTTTTCTTCCAGCACCAGCAGCCCGTAACT
>a8d321b9-97f1-42aa-a17a-fa1c4609e60e runid=b48dcb2029ef589b8cdb55fe529d93b3baadbc7c sampleid=28012021COVID3 read=46034 ch=83 start_time=2021-01-28T18:24:16Z
CGTTCGGTGCGTATTGCTGGTGCTGAAGAAAGTTGTCCGGTGTCTTTGTGTTAACCTTTCTGTTGGTGCTGATATTGCGGCGTCTGGCAGGTGTTAACCTGGCCTCGAGAGAGTTTGATCCTGGCTGGGATGAACGCTGACGTGCCTAATACATGCAAGTAG
>6d97cca8-6795-4c59-b4d1-05880735b957 runid=b48dcb2029ef589b8cdb55fe529d93b3baadbc7c sampleid=28012021COVID3 read=46363 ch=227 start_time=2021-01-28T18:24:08Z
GTTCAGTGCATATTGCTGTGGTGCTGCGAAAGTTTGTGGTGTCTTTGTGTTAACGCTTTCTGTTGGTGCTGATATTGCGGCGGCGTCTGCGCGGTGTTTAACCTGGCGGATCAAGGAGGTGTTCAGCCGCACCTTCCGTTTTGACTACGCCATTACGACTGCCCAATCTGTTTTACCCTAGGCCCGATCCTTGCGGTCACGGACTTCAGGCACCCCCGGCTTTCATGGCTTGACAGGCGGTGTGTACAAGGCCCGGGAACATTCACCCGCGCCATGGCTGATGCGCGATTACTAGCGAATCCCAGCTTCACGGGTCGGGTTGCA
>772e0bef-fad4-4767-b0e6-d88e6f16b2ed runid=b48dcb2029ef589b8cdb55fe529d93b3baadbc7c sampleid=28012021COVID3 read=49763 ch=349 start_time=2021-01-28T18:24:10Z
CGTTCAGTTACGTATTGCTGGTGCTGAAGAAAGTTGTCGGTGTCTTTAGGTTAACCTTTCTGTTGGTGCTGATATTGCGGCGTCTGCTTGGGTGTTTAACCTGGCCTCGAGAAGTTTGATCCTGGCTCAGGATGAACGCTAGCTTCTGAGCTTAACATGCAGAGTCCAGGGCAGCATGGAAGAAACTGCTTCTTCTGATGGCGACAACGCACGGGTGCATGCGCGTATCAAACCTGCCTCATA
>2b673c48-69c9-48a4-995b-4ad630e23502 runid=b48dcb2029ef589b8cdb55fe529d93b3baadbc7c sampleid=28012021COVID3 read=71966 ch=124 start_time=2021-01-28T18:24:11Z
TTCGTTCAGTTACGTATTGCTGGTGCTGAAAGTTGTCGGTGTCTTTGTGTTAACCTTTCTGTTGGTGCTGATATTGCGGCTCTGCTTGGGTGTTTAACCTGGCGGATCAAGGAGTGTTAGCCGCACCTTCCGGTACCAACTACCTTGTTACGACTTAGCCCC

OTUs table:
OTU ID  IA1 IA2 IA3 IA4 IA5 IS6
otu500  7   10  3   0   5   0
otu502  7   8   4   0   11  1
otu504  4   5   6   2   4   2
otu505  4   4   6   2   8   3

and my desired output is that the .fasta file would look like this:
>otu500
tggggaatcttagacaatgggcgcaagcctgatctagccatgccg
>otu502
tgaggaatattggtcaatggaggcaactctgaaccagccatgccttggtcaatggaggc
>otu504
gcgaacaggattagataccctggtagtccacgccgtaaacttggtcaatggaggc
>otu505
ttggtcaatggaggcttggtcaatggaggctaccctggttaccctggt
>otu506

There is not enough information here to label the fasta sequences with their OTU assignments. What program did you use to produce the OTU table? The program should have a way to list each sequence's OTU assignment. For example, in the mothur program, the cluster command outputs a list file.

There seems to be some missing info - I can't see any logical mapping between reads and OTU table - do you have an OTU sequence file?


text:
Naming functions in R
I have a set of files based on which I expect to load RDS files like R1-Tn5-NHP-55_NA and R1-Tn5-NHP-42_NA but the code is producing just the name of the 55 files not both 55 and 42.
#standardize N set names
rename_KP_set <- function(sample_name_vector)  {
 new_sample_name_vector <- apply(sample_name_vector,
      function(x){
         if (grepl("42",x[3])){
           c(x[1:3],'55',x[5])
         } else if (grepl("55",x[3])) {
           c(x[1:3],'NHP',x[5])
         } else if (grepl("55",x[3])){
           c(x[1:3],'NHP',x[5])
 } else
c(x[1:3],'55',x[5])
 }
}) %>%
     sapply(paste, collapse='-')
     return(new_sample_name_vector)
 }

How I can get separate rds files having 42 and 55 numbers or for that matter any rds with some number attached to them?
I have already posted this question on Stack overflow on February 27th, 2023 but have got no answers so I am reposting that in Stack exchange bioinformatics in search of an answer. Any help on this issue will be considered useful.

Your question has a few wording issues. For example, what do you mean by "...a set of files based on which I expect to load RDS files..."? Provide sample input, expected output and actual output. You'll probably figure out the problem on your own as you prepare these items.

Please [edit] your question and provide us with input files and expected output so we can reproduce the problem and test our solutions. We cannot work on a problem if we cannot reproduce it.

Can you please explain what your function is trying to do? It's not clear to me from just looking at the function.


text:
Cut and Tag analysis
I am trying to do Cut and Tag analysis for a set of files and I want to perform a procedure called Tn5 normalization of samples.
The basic goal is to generate bedgraph files and does normalization with the weights of a sequencing depth. I have a set of function files that I have defined and then a set of Tn5 normalization scripts. The scripts are both quite long and hence I will be sharing parts of the scripts to get hold of my error and also provide as much information as possible.
My bedgraph files look like this:
replicate species Msample serotype protein            sample      condition    group
1        R1     NHP      02     LK03      NA R1-NHP-02-LK03-NA R1-NHP-LK03-NA LK03-NHP
2        R1     NHP      04     LK03      NA R1-NHP-04-LK03-NA R1-NHP-LK03-NA LK03-NHP
3        R1     NHP      13      KP1      NA R1-NHP-13-KP1-NA R1-NHP-KP1-NA KP1-NHP
4        R1     NHP      42      KP1      NA R1-NHP-42-KP1-NA R1-NHP-KP1-NA KP1-NHP
5        R1     NHP      55     LK03      NA R1-NHP-55-LK03-NA R1-NHP-LK03-NA LK03-NHP
6        R1     NHP      57      KP1      NA R1-NHP-57-KP1-NA R1-NHP-KP1-NA KP1-NHP

I must load a complete set of files called RDS which basically bears a list of the individual bed graph files.
The read lengths for each of the files look like this:
seqnames start  end normcount
73694699     opt2     0   52 0.0000000
73694700     opt2    52   56 0.0243838
73694701     opt2    56   83 0.0731515

The Tn5 R normalization scripts which are relevant
library(dplyr, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")
library(gridExtra, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")
library(Gviz, lib.loc = "/home/ag111/R/x86_64-pc-linux-gnu-library/4.1") #needed to be installed
library(ggplot2)
library(RcppArmadillo, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")
#library(stringi, lib.loc = "/scg/apps/software/r/4.0.3/lib64/R/library")
library(tidyr, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")
library(RColorBrewer, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")
library(ggpubr, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library") #needed to be installed
library(GenomicFeatures, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")
library(GenomicRanges, lib.loc = "/scg/apps/software/r/4.1.2/lib64/R/library")

source('AG_CnRT_functions.r')
dir.create(file.path("Tn5exITR_InVivo_norm_full"))

intermediates_dirs <- c(SeqDepth='intermediates/')
#intermediates_dirs <- c(SeqDepth='/labs/markay/Aranyak/Tn5_NHP_rename/Data_analysis/intermediates/')

normalization_method <- 'SeqDepth'

AAV_gsize <- 4410
AAV_anno <- read.delim('AAV_opt2_anno.txt')

#figures_dir <- c(SeqDepth='Tn5exITR_InVivo_norm_full2/')

figures_dir <- c(SeqDepth='Tn5exITR_InVivo_norm_full2/')

############## Tn5 ###########
AAV_normbedgraphs <- readRDS("/labs/markay/Aranyak/Tn5_NHP_rename/Data_analysis/intermediates/AAV_rmdup_SeqDepthNormbedgraphs.rds")

naming_order <- c('replicate','species',"Msample",'serotype','protein')

AAV_normbedgraphs_exITR <- AAV_normbedgraphs %% lapply(function(x){
  df <- x %% subset(start = 185 & end <=  4167)
  return(df)
})

AAV_normcov_regions_meta_exITR <- list()
AAV_normcov_regions_meta_exITR[[normalization_method]] <- get_normcov_regions_meta_AG(
  AAV_normbedgraphs=AAV_normbedgraphs_exITR,
  naming_order=naming_order,
  AAV_gsize=AAV_gsize,
  AAV_anno=AAV_anno)

Merging data
AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype[AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype == "LK"] <- "LK03"
AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype[AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype == "iG"] <- "AM"
AAV_normcov_regions_meta_exITR[naming_order] <- lapply(AAV_normcov_regions_meta_exITR[naming_order], factor)
AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype <- ordered(AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype, levels = c("AM", "LK03", "DJ"))

AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype[AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype == "LK"] <- "LK03"
AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype[AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype == "iG"] <- "KP1"
AAV_normcov_regions_meta_exITR[naming_order] <- lapply(AAV_normcov_regions_meta_exITR[naming_order], factor)
AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype <- ordered(AAV_normcov_regions_meta_exITR[[normalization_method]]$serotype, levels = c("KP1", "LK03"))

metadata <- get_metadata_AG3(AAV_normbedgraphs=AAV_normbedgraphs_exITR,
                             #naming_order=naming_order)
metadata <- get_metadata_KP(AAV_normbedgraphs=AAV_normbedgraphs_exITR,
                             naming_order=naming_order)

head(metadata)

#coverage separated by AAV regions and species
#to_plot <- c()
#to_plot <- AAV_normcov_regions_meta_exITR[[normalization_method]] %% gather('region','normcov',2:6)
#to_plot$region <- factor(to_plot$region, levels=c(as.character(AAV_anno_exITR$region),'AAV_normcov'))
#max_y <- ceiling(max( log2(to_plot$normcov1) ) )
#min_y <- floor( min( log2(to_plot$normcov1) ) ).

After loading of the rds file the error I am getting while I am executing the code block
 AAV_normcov_regions_meta_exITR <- list()
 #AAV_normcov_regions_meta_exITR[[normalization_method]] <- get_normcov_regions_meta_AG3(
 AAV_normcov_regions_meta_exITR[[normalization_method]] <- get_normcov_regions_meta_AG(
   AAV_normbedgraphs=AAV_normbedgraphs_exITR,
   naming_order=naming_order,
   AAV_gsize=AAV_gsize,
   AAV_anno=AAV_anno)

Error in if (!all(AAV_normcov_regions$sample == metadata$sample)) { :
  missing value where TRUE/FALSE needed.

The code block which I think for the function file that is valid  is
AAV_normcov_regions[is.na(AAV_normcov_regions)] <- 0
  AAV_normcov_regions[is.null(AAV_normcov_regions)] <- 0
  AAV_normcov_regions[is.na(AAV_normcov_regions)] = 0
  AAV_normcov_regions <- AAV_normcov_regions[ match(metadata$sample,AAV_normcov_regions$sample) , ]
  if ( !all(AAV_normcov_regions$sample == metadata$sample) ){
    stop('Unsuccessful matching to metadata')
  }
  return(AAV_normcov_regions).

This function is working for another set of samples to generate a plot list.
The plot list I am having is to_plot
[1] sample      AAV_normcov replicate   species     Msample     serotype    protein     condition   group
<0 rows>

(or 0-length row.names) where as what is expected is something like this
  sample   Enhancer        FLuc      polyA   Promoter       WPRE AAV_normcov replicate batch species drug serotype condition group
1  R1-LiD-H-N-DJ 0.42244801 0.290929537 0.43640203 0.22379898 0.34125606  0.26650388        R1   LiD       H    N       DJ      H-DJ  H-DJ
2  R1-LiD-H-N-iG 1.42930414 0.259549496 0.57248120 0.37425687 0.42189340  0.36092640        R1   LiD       H    N       AM      H-iG  H-iG
3  R1-LiD-H-N-LK 0.25754367 0.317786560 0.27373645 0.19470267 0.29818285  0.24753280        R1   LiD       H    N     LK03      H-LK  H-LK
4  R1-LiD-M-N-DJ 0.24578258 0.097169668 0.14489638 0.06470209 0.11091813  0.09308374        R1   LiD       M    N       DJ      M-DJ  M-DJ

I am happy to share other code blocks if it is required to answer the question.

Please see our help page on formatting. I tried to format your code as code, but please check and make sure I didn't miss something.

Is there a typo in your code? You have a x %% subset which should probably be x %>% subset. Did the code run or did you edit it after copy-pasting and make a mistake there?

It will be indeed be x %>% subset. There was a mistake during copy pasting.

Check your code for more copy-paste errors and make sure everything lines up. Edit your post and fix all the errors you find, including the %% subset error. Please use proper code formatting for inline code and code blocks.

The file that you include as a bedgraph example doesn't look like bedgraphs.


text:
Adapter trimming
I am trying to do adapter trimming, alignment and sorting for a range of large scale  paired end fastq files. The code I am using is given below:
#!/bin/bash
# Do Trimming, Alignment according to Config file

#SBATCH --job-name=T_Align_single
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=32
#SBATCH --partition=batch
#SBATCH --account=markay
#SBATCH --time=72:00:00
#SBATCH --mem=16G
#SBATCH --mail-type=TIME_LIMIT_80,BEGIN,END
#SBATCH --mail-user=<email-redacted>

module load bowtie2/2.4.1
module load samtools/1.9

mkdir -p Mapping/

config=${INPUT_FILE}

file1=$(cat $config | awk '$1 ~ /fastq/ {print $1}')
file2=$(cat $config | awk '$1 ~ /fastq/ {print $2}')
file3=$(cat $config | awk '$1 ~ /index=/ {print $1}'| sed 's/ref_index=//')

#adapter trimming
trimfile1=Mapping/"${file1%.fastq}".trim.fastq
trimfile2=Mapping/"${file2%.fastq}".trim.fastq

 java -jar /labs/markay/Aranyak/Adriana_intro_CnT/Cnt_Test/Test-analysis/fastq_processing/Trimmomatic-0.39/trimmomatic-0.39.jar  PE -phred33 $file1 $file2 \
   $trimfile1 Mapping/"${file1%.fastq}".unpaired.fastq \
   $trimfile2 Mapping/"${file2%.fastq}".unpaired.fastq \
   LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 \
   MINLEN:20 ILLUMINACLIP: /labs/markay/Aranyak/Adriana_intro_CnT/adapter.fa:2:30:10 -threads ${SLURM_CPUS_PER_TASK}

# mapping using bowtie2
sortbam=Mapping/"${file1%_R1_001.fastq}".sort.bam
sample_ID=$(cat $config | awk '$1 ~ /fastq/ {print $4}')
bowtie2 -p ${SLURM_CPUS_PER_TASK} --local --very-sensitive-local --no-unal \
   --no-mixed --no-discordant --phred33 -I 10 -X 700 \
   -x $file3 \
   -1 $trimfile1 -2 $trimfile2 --rg-id $sample_ID | samtools view -Sb - | samtools sort -@${SLURM_CPUS_PER_TASK} -o $sortbam

#Index the sorting by coordinate bam file
samtools index ${sortbam}

My original fastq files before adapter trimming looked like this
R3-LiD-M-LK_S16_L002_R2_001.fastq.gz
R3-LiD-M-LK_S16_L002_R1_001.fastq.gz

However after adapter trimming I removed the lane names and my file looks like this
R3-LiD-M-LK_S16_R2_001.fastq.gz
R3-LiD-M-LK_S16_R1_001.fastq.gz

The code that worked for the original file is not generating a proper bam file at the end though I am getting a list of trimmed fastq files. The error message I am getting is:
[E::hts_open_format] Failed to open file R1-LiD-H-DJ_S3.sort.rmdup.bam
samtools index: failed to open "R1-LiD-H-DJ_S3.sort.rmdup.bam": No such file or directory
[E::hts_open_format] Failed to open file R1-LiD-H-DJ_S3.sort.bam
The file 'R1-LiD-H-DJ_S3.sort.bam' does not exist
[E::hts_open_format] Failed to open file R1-LiD-H-DJ_S3.sort.rmdup.bam
The file 'R1-LiD-H-DJ_S3.sort.rmdup.bam' does not exist
R1-LiD-H-DJ_S3 DONE.

How should I change the code to generate proper alignment and bam files?

This should be a simple learning exercise for you - all that's changed is that the lane number part is gone, but your error doesn't seem related to that at all, and no variables use the lane number part anywhere. Add a set -x after the #SBATCH lines and look at the job's STDERR file to see the exact commands being executed.

Or just print out your variables so you can see what is happening. You are trying to open a file named R1-LiD-H-DJ_S3.sort.bam and a file named R1-LiD-H-DJ_S3.sort.rmdup.bam and neither of them exist. You also don't mention rmdub anywhere in your code, is this really the code you are running?

Thanks for your help, I have been able to trouble shoot the error with the naming of the files.

Hi @AranyakGoswami, can you please write up your solution as an answer to this question (and make it an accepted answer)? That will help the StackExchange bots do the right thing with regards to not hunting for more answers to an already-solved question.


text:
Parallelizing microRNA targets
I  am trying to look for miRNA targets using a file called Conservedfamily.txt from the Zebrafish target scan fish website. I have written a python program to extract each of the miRNAs which works fine. However, for a big list of microRNAs, I must provide a name every time in line 22 and get one individual output file. I am thinking is there any option to parallelize the process so that somehow, I can provide a list of all microRNAs and obtain all my targets at the same time. The python code is attached below. I have posted the question in biostars but go no reply from them for 4 days. Hence I am posting the question in Stack exchange. The python code is attached below.
#!/usr/bin/python

#Import relevant modules

import sys

conservedfamily = open(sys.argv[1], 'r')
targetgenelist = open(sys.argv[2], 'w')
targetgeneid = open(sys.argv[3], 'w')

geneid = []

targetgenelist.write('miR family' + '\t' + 'Gene ID' + '\t' + 'Gene Symbol' + '\t' + 'Transcript ID' + '\t' + 'Species ID' + '\t' + 'UTR start' + '\t' + 'UTR end' + '\t' + 'MSA start' + '\t' + 'MSA end' + '\t' + 'Seed match' + '\n')

for line in conservedfamily:
  splitline = line.strip().split('\t')

  if splitline[0] == 'miR-210':
    targetgenelist.write(line)
    if not splitline[1] in geneid:
      geneid.append(splitline[1])
      targetgeneid.write(splitline[1] + '\n')

conservedfamily.close()
targetgenelist.close()
targetgeneid.close()

python miR99target.py 'Conserved_Family_Info.txt' 'miR-210targetlist.txt' 'miR-210targetid.txt'

I'd suggest loading in all of the miRNA tags as a set(). Then, you can just check the set membership each time.
#!/usr/bin/python

#Import relevant modules

import sys

conservedfamily = open(sys.argv[1], 'r')
targetgenelist = open(sys.argv[2], 'w')
targetgeneid = open(sys.argv[3], 'w')

# Create a set of all miRNA tags. Each tag should be on its own line.
with open(sys.argv[4], 'r') as miRNAlist:
    miRNA_set = {line.strip() for line in miRNAlist}

geneid = []

targetgenelist.write('miR family' + '\t' + 'Gene ID' + '\t' + 'Gene Symbol' + '\t' + 'Transcript ID' + '\t' + 'Species ID' + '\t' + 'UTR start' + '\t' + 'UTR end' + '\t' + 'MSA start' + '\t' + 'MSA end' + '\t' + 'Seed match' + '\n')

for line in conservedfamily:
  splitline = line.strip().split('\t')

  # Check if current tag is in the desired set
  if splitline[0] in miRNA_set:
    targetgenelist.write(line)
    if not splitline[1] in geneid:
      geneid.append(splitline[1])
      targetgeneid.write(splitline[1] + '\n')

conservedfamily.close()
targetgenelist.close()
targetgeneid.close()

One additional comment:
I'd suggest using the python module for csv reading/writing, even for quick and dirty scripts as well as using context managers for reading/writing files.


text:
Setting Contrast for DESeq2 results
This question was also asked on BioStars
I am trying to recreate this heatmap. How would I compare the four variables: Ly49+/- and MOG/MOGSP to get a single heatmap? Thank you for helping me through this.
I need to compare both Ly49+ versus Ly49− and MOG versus MOG plus SP. Column data in at the bottom. The function I used for DESeq2 is,
dds = DESeqDataSetFromMatrix(countData=countData,
                                colData=colData,
                                design= ~ cell_treatment)

I am not sure if the design is right. For the results, how should I write the contrasts? I am very confused and would really appreciate any help!
DataFrame with 11 rows and 2 columns
cell treatment
<character> <character>
MOGSP_2 SPtetramer+CD8+TCell MOGSP
MOGSP_3 SPtetramer+CD8+TCell MOGSP
MOG_3 SPtetramer+CD8+TCell MOG
MOG_4 SPtetramer+CD8+TCell MOG
MOG_5 SPtetramer+CD8+TCell MOG
Ly49N_1 Ly49-CD8+TCell MOGSP
Ly49P_1 Ly49+CD8+TCell MOGSP
Ly49N_2 Ly49-CD8+TCell MOGSP
Ly49P_2 Ly49+CD8+TCell MOGSP
Ly49N_3 Ly49-CD8+TCell MOGSP
Ly49P_3 Ly49+CD8+TCell MOGSP

Thank you!

Your table is not very clear. If I understood it correctly, there are three groups of t-cells, SPtetramer+, Ly49- and Ly49+ and within it there is a MOGSP vs MOG comparison. I would set colData like this:
                   name             Ly treatment
1  SPtetramer+CD8+TCell SPtetramerplus     MOGSP
2  SPtetramer+CD8+TCell SPtetramerplus     MOGSP
3  SPtetramer+CD8+TCell SPtetramerplus       MOG
4  SPtetramer+CD8+TCell SPtetramerplus       MOG
5  SPtetramer+CD8+TCell SPtetramerplus       MOG
6        Ly49-CD8+TCell      Ly49minus     MOGSP
7        Ly49+CD8+TCell       Ly49plus     MOGSP
8        Ly49-CD8+TCell      Ly49minus     MOGSP
9        Ly49+CD8+TCell       Ly49plus     MOGSP
10       Ly49-CD8+TCell      Ly49minus     MOGSP
11       Ly49+CD8+TCell       Ly49plus     MOGSP

I generate an example dataset:
counts = matrix(rnbinom(11000,mu=100,size=1),ncol=11)

And in the design you specify these two factors, meaning gene expression can be affected by these two conditions:
dds = DESeqDataSetFromMatrix(counts,colData=colData,~Ly +treatment)
dds = DESeq(dds)

You can compare them like this:
head(results(dds,c("Ly","Ly49plus","Ly49minus")))
log2 fold change (MLE): Ly Ly49plus vs Ly49minus
Wald test p-value: Ly Ly49plus vs Ly49minus
DataFrame with 6 rows and 6 columns
          baseMean     log2FoldChange            lfcSE               stat
         <numeric>          <numeric>        <numeric>          <numeric>
1 91.4458114779915  0.712987165797688 1.28291188260395  0.555756927241586
2 81.7945234843954 -0.418266530066077 1.14378290088875 -0.365686993345568
3 94.4351681904941 -0.775592496878288 1.32965355129901 -0.583304196886903
4 107.484335203298  -1.36694885973712 1.38595911743023 -0.986283680770933
5 61.2205696982841 -0.760219864564472 1.16779888688769 -0.650985262188887
6 95.5908168437137  0.497427543096335 1.12768919469662   0.44110340458672
             pvalue              padj
          <numeric>         <numeric>
1 0.578377034355333 0.977537706056193
2    0.714598652706 0.992321035598301
3 0.559688537757548 0.977537706056193
4 0.323993925580938 0.954772990635094
5 0.515056000552832 0.977537706056193
6 0.659138138867294 0.992321035598301

To compare groups within another factor:
results(dds,c("treatment","MOGSP","MOG"))
log2 fold change (MLE): treatment MOGSP vs MOG
Wald test p-value: treatment MOGSP vs MOG
DataFrame with 6 rows and 6 columns
          baseMean     log2FoldChange            lfcSE               stat
         <numeric>          <numeric>        <numeric>          <numeric>
1 91.4458114779915 -0.875905107403943 1.43007567177694 -0.612488642867122
2 81.7945234843954  0.595588650280833 1.28109080878515  0.464907441530728
3 94.4351681904941 0.0255578400583986 1.48592401497315 0.0171999643325371
4 107.484335203298  -1.46652057431159 1.55285829334496 -0.944400774105796
5 61.2205696982841  -1.27442596494316 1.30489822606005  -0.97664778715433
6 95.5908168437137  -2.35402796630456 1.26608624796739  -1.85929510733078
              pvalue              padj
           <numeric>         <numeric>
1  0.540214509839105 0.944881849859735
2  0.641997741235223 0.973651054534338
3  0.986277090644156 0.994944323499595
4  0.344964886056295  0.87421063696937
5  0.328743552427222                NA
6 0.0629853201285643 0.525889480109579

The colData i used:
structure(list(name = structure(c(3L, 3L, 3L, 3L, 3L, 1L, 2L,
1L, 2L, 1L, 2L), .Label = c("Ly49-CD8+TCell", "Ly49+CD8+TCell",
"SPtetramer+CD8+TCell"), class = "factor"), Ly = c("SPtetramerplus",
"SPtetramerplus", "SPtetramerplus", "SPtetramerplus", "SPtetramerplus",
"Ly49minus", "Ly49plus", "Ly49minus", "Ly49plus", "Ly49minus",
"Ly49plus"), treatment = structure(c(2L, 2L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L), .Label = c("MOG", "MOGSP"), class = "factor")), class = "data.frame", row.names = c(NA,
-11L))

so it can used for multiple variable or cell type comparisons instead of many pairwise comparison which i do normally

it depends on the design but in general yes. In this example, you can see it's necessary because the groups have some overlap

conceptual doubt so if i have to perform multiple pairwise comparison do i need to set reference level for every comparison or i can simply use "results(dds,c("treatment","MOGSP","MOG"))" or the contrast function?

using results(dds,c("treatment","MOGSP","MOG")) takes care of everything. The reference level of the factor does not matter here

will try this..


text:
Fetching Protein Sequences through UniPROT gives ConnectionError
I have 193000 protein interactions in CSV named all_proteininteractions.csv which have query protein name in the 2nd column and partner protein in the 3rd like this:
Query_ENSP,Query_Name,Partner_Name,Combined_Score
9606.ENSP00000000233,ARF5,RAC2,0.4
9606.ENSP00000000233,ARF5,RAC2,0.4
9606.ENSP00000000233,ARF5,PPP2CB,0.408
9606.ENSP00000000233,ARF5,PPP2CB,0.408
9606.ENSP00000000233,ARF5,R3HDML,0.413
9606.ENSP00000000233,ARF5,R3HDML,0.413
9606.ENSP00000000233,ARF5,PLEK,0.422
9606.ENSP00000000233,ARF5,PLEK,0.422
9606.ENSP00000000233,ARF5,RAB11FIP5,0.501
9606.ENSP00000000233,ARF5,RAB11FIP5,0.501
9606.ENSP00000000233,ARF5,AP3B1,0.553
9606.ENSP00000000233,ARF5,AP3B1,0.553
9606.ENSP00000000233,ARF5,CYTH4,0.579
9606.ENSP00000000233,ARF5,CYTH4,0.579
9606.ENSP00000000233,ARF5,RABAC1,0.658
9606.ENSP00000000233,ARF5,RABAC1,0.658
9606.ENSP00000000233,ARF5,YKT6,0.672
9606.ENSP00000000233,ARF5,YKT6,0.672
9606.ENSP00000000233,ARF5,BET1,0.679

I am using this Python code to fetch sequences of each protein and store them in a CSV file:
import pandas as pd
import requests
import math

requests.packages.urllib3.disable_warnings()

# Function to fetch protein sequence from UniProt
def fetch_protein_sequence_from_uniprot(protein_name):
    #Declaring UniProt API
    uniport_api_url = f"https://rest.uniprot.org/uniprotkb/search?query=gene_exact:{protein_name}+AND+organism_id:9606&format=fasta&size=1"
    response = requests.get(uniport_api_url,verify=False)

    # Parse the response to extract sequence
    sequence = ""
    if response.ok:
        lines = response.text.split("\n")
        for line in lines:
            if not line.startswith(">"):  # Exclude header lines
                sequence += line
    return sequence

# Reading the input CSV file
input_csv = "all_protein_interactions.csv"
chunksize = 1000  # Batch size

# Create an iterator to read the CSV in chunks
csv_reader = pd.read_csv(input_csv, chunksize=chunksize)

# Initialize an empty list to store modified chunks
modified_chunks = []

# Process each chunk
for chunk in csv_reader:
    # Fetch protein sequences for query and partner names
    query_sequences = [fetch_protein_sequence_from_uniprot(name) for name in chunk.iloc[:, 1]]
    partner_sequences = [fetch_protein_sequence_from_uniprot(name) for name in chunk.iloc[:, 2]]

    # Add fetched sequences as new columns in the chunk
    chunk['Query_Sequence'] = query_sequences
    chunk['Partner_Sequence'] = partner_sequences

    # Append the modified chunk to the list
    modified_chunks.append(chunk)

# Concatenate all modified chunks into a single DataFrame
result_df = pd.concat(modified_chunks, ignore_index=True)

# Save the updated DataFrame back to CSV
result_df.to_csv("protein_interactionswithseqs.csv", index=False)

First, it was throwing SSLError, so I made verify=False. But now, after running the script for 3 hours it's showing following error:
ConnectionError: ('Connection aborted.', RemoteDisconnected('Remote end closed connection without response'))

How can I remove this connection aborted error to fetch 386000 sequences?

You also should add print statements in various places, to see where the problem is. Has the script managed to finish the for chunk loop? Where exactly is it failing? And isn't your loop trying to get the same protein multiple times? You seem to have the same name in the 2nd column for several lines, so why get that same protein multiple times?

Thanks. I am using  Windows 10 with Intel(R) Core(TM) i7-3770 CPU and 12 GB ram.

I modified the text as well. I have similar proteins but I want to fetch them all and store in 5th and 6th column

Your script works as expected on the example you gave. Running your script results in a protein_interactionswithseqs.csv with lines like:
$ head -n2 protein_interactionswithseqs.csv
Query_ENSP,Query_Name,Partner_Name,Combined_Score,Query_Sequence,Partner_Sequence
9606.ENSP00000000233,ARF5,RAC2,0.4,MGLTVSALFSRIFGKKQMRILMVGLDAAGKTTILYKLKLGEIVTTIPTIGFNVETVEYKNICFTVWDVGGQDKIRPLWRHYFQNTQGLIFVVDSNDRERVQESADELQKMLQEDELRDAVLLVFANKQDMPNAMPVSELTDKLGLQHLRSRTWYVQATCATQGTGLYDGLDWLSHELSKR,MQAIKCVVVGDGAVGKTCLLISYTTNAFPGEYIPTVPGRPCALLQDAEAGLSSDVEGACWCQTQEGRMTGEAPAPPWPALRGPHHS

However, this is horribly, horribly inefficient! You are fetching the same sequence multiple times. With your example data, the protein ARF5 is actually downloaded 19 times, meaning this will take 18 times more time than needed. This is very likely the cause for the timeout error you describe. Now, it might be that your sequences are too many to be stored in RAM (but I would still try it), so one way around this is to download the sequences first, and then write the mappings. Something like this:
import pandas as pd
import requests
import sys
requests.packages.urllib3.disable_warnings()

def write_protein_seq_file(protein_name, protein_seq):
    file_name = "seqs/" + protein_name
    seq_file = open(file_name, "w")
    data = f"{protein_seq}\n"
    seq_file.write(data)
    seq_file.close

def read_protein_seq_file(protein_name):
    seq_file = open("seqs/" + protein_name, "r")
    data = seq_file.read()
    seq_file.close
    return data.rstrip()

# Function to fetch protein sequence from UniProt
def fetch_protein_sequence_from_uniprot(protein_name):
    #Declaring UniProt API
    uniport_api_url = f"https://rest.uniprot.org/uniprotkb/search?query=gene_exact:{protein_name}+AND+organism_id:9606&format=fasta&size=1"
    response = requests.get(uniport_api_url,verify=False)

    # Parse the response to extract sequence
    sequence = ""
    if response.ok:
        lines = response.text.split("\n")
        for line in lines:
            if not line.startswith(">"):  # Exclude header lines
                sequence += line
    return sequence

# Reading the input CSV file
input_csv = "all_protein_interactions.csv"
chunksize = 1000  # Batch size

# Create an iterator to read the CSV in chunks
csv_reader = pd.read_csv(input_csv, chunksize=chunksize)

# Initialize an empty list to store modified chunks
modified_chunks = []

# Get all sequences
for chunk in csv_reader:
    for protein in set(list(chunk.iloc[:, 1]) + list(chunk.iloc[:, 2])):
        protein_seq = fetch_protein_sequence_from_uniprot(protein)
        write_protein_seq_file(protein,protein_seq)

print("Finished fetching sequences", file=sys.stderr)
## Now, process the file again and write the new file
csv_reader = pd.read_csv(input_csv, chunksize=chunksize)

for chunk in csv_reader:

    # Fetch protein sequences for query and partner names
    query_sequences = [read_protein_seq_file(name) for name in chunk.iloc[:, 1]]
    partner_sequences = [read_protein_seq_file(name) for name in chunk.iloc[:, 2]]

    # Add fetched sequences as new columns in the chunk
    chunk['Query_Sequence'] = query_sequences
    chunk['Partner_Sequence'] = partner_sequences

    # Append the modified chunk to the list
    modified_chunks.append(chunk)

# Concatenate all modified chunks into a single DataFrame
result_df = pd.concat(modified_chunks, ignore_index=True)

# Save the updated DataFrame back to CSV
result_df.to_csv("protein_interactionswithseqs.csv", index=False)

First, create a directory named seqs, then run the script above. It will process the csv file twice. Once to get all sequences mentioned in the file and save each one in a file called seqs/$sequence_name (where $sequence_name is replaced by the actual name of each protein. Then, it goes through the file a second time, but instead of the very time consuming download, it now fetches the sequence from the file stored on your computer. Once you are satisfied this works, you can delete the seqs directory.
Note that I don't have access to a Windows machine, so this was tested on Linux. Also, I am not a python coder, so I am sure this ins't very good python and there will be improvements to be made, but even on the toy example you gave, this speeds things up significantly:
$ time your_script.py

real    0m10.186s
user    0m0.998s
sys     0m0.076s

$ time my_script.py
Finished fetching sequences

real    0m3.180s
user    0m0.456s
sys     0m0.066s

As you can see above, even on the tiny example with only 19 data lines, your script took 10 seconds while mine only took 3. So try this approach (or, better, take my example and make it more pythonic, add some error checking, some try/except blocks, make sure your data are correct) and you might find the timeout error is no longer an issue. You should certainly find that the whole thing takes far less than three hours(!) to run.
First of all, that is a really, really inefficient way of doing it. Using your example data, you are getting the same protein, ARF5, 19 times. This makes your program 18 times slower than it needs to be. You should modify the script to download each protein once and then write them all out.

Your example script is not taking any file and where to download the all sequences

@AnasJamshed what do you mean? It is taking the exact same file as your version and in the exact same way. The sequences will be downloaded into the seqs directory which, as I explain in my answer, you need to create.


text:
How to Split VCF File into male and Female groups
Cross-posted on biostars
I have trio_example.vcf file and I want to split them into male and female groups. I tried this:
$ bcftools query -l trio_example.vcf

And it gave me:
N1
N2
N3

After that I tried this command:
for sample in `bcftools query -l trio_example.vcf.gz`; do vcf-subset --exclude-ref -c $sample trio_example.vcf.gz > ${sample}.vcf; done

to split 3 samples. Now its showing same data in all 3 N files
How can I solve this issue? Also I need to split the data into male and female. For splitting I am using:
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%AF\n' trio_example.vcf | awk '{if($1 == "X" && $5 >= 0.3 && $5 <= 0.7) {print $6 " is male"} else if ($1 == "Y" && $5 == 1) {print $6 " is male"} else if ($1 == "X" && $5 == 1) {print $6 " is female"} else {print $6 " is of unknown sex"}}'

But it also not give me exact result

Do you know which of the individuals in the VCF are male and female? If so, you can use bcftools (view command) to subset the VCF twice, once creating a VCF only including males, a second time making a VCF only including females.

I dont know thats why I am using now this command: bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%AF\n' trio_example.vcf | awk '{if($1 == "X" && $5 >= 0.3 && $5 <= 0.7) {print $6 " is male"} else if ($1 == "Y" && $5 == 1) {print $6 " is male"} else if ($1 == "X" && $5 == 1) {print $6 " is female"} else {print $6 " is of unknown sex"}}'

Is the information about the sex of your individuals anywhere in the VCF? In the INFO tag maybe? Somewhere in any paper or supplementary file? Which organism are you working on -- does this organism have sex chromosomes? Does the VCF have X and Y chromosomes (or Z/W etc depending on the sex determining system)? If so, individuals with variants at Y or hemizygous for X should be male, whereas individuals with homozygous variants at X should be female.

I only have X chromosome present in VCF File

What's the species. Drosophila?

I dont know about the species but when I separate N1 , N2 and N3 all files show the same record

If you don't even know what species you're working with I think you should figure that out ... because if you want to make an inference of sex from the VCF and that information simply isn't stated anywhere, you need to know the sex determining system you're working with. Drosophila only has X chromosomes, but males have one X and females have two. This kind of information is important if you want to infer sex because males can only have one variant at X, females can have two.

@AnasJamshed please add all extra information to your question, not in comments. As you can see your bcftools -f command is  very hard to read. Add it to the question please, along with all the other extra information you have given so that the next person to see the question can have all information there.

@arara don't forget the pseudoautosomal regions: you will likely have some, albeit few, homozygous variants in both X and Y even in males.

Cross-posted on biostars. OP was reprimanded on biostars by multiple users for opening a lot of tangential posts and is probably attempting a similar move here.


text:
Making a box plot of the proportion of cells in each cell types in two groups
I have number of cells in three cell types T, B, M in a Seurat object
For two groups of patients, cancers and controls like
> table(Idents(pbmc), pbmc$samples)

    PN0252_0001 PN0252_0002 PN0252_0003 PN0252_0004 PN0252_0005 PN0252_0006 PN0252_0008
  T        3550         568         136         268        1835           2        1678
  B         136         129          11          25         133           0         203
  M         120          22          18          26         159           0          59
> dput(table(Idents(pbmc), pbmc$samples))
structure(c(3550L, 136L, 120L, 568L, 129L, 22L, 136L, 11L, 18L,
268L, 25L, 26L, 1835L, 133L, 159L, 2L, 0L, 0L, 1678L, 203L, 59L
), .Dim = c(3L, 7L), .Dimnames = structure(list(c("T", "B", "M"
), c("PN0252_0001", "PN0252_0002", "PN0252_0003", "PN0252_0004",
"PN0252_0005", "PN0252_0006", "PN0252_0008")), .Names = c("",
"")), class = "table")
>

Patients PN0252_0001, PN0252_0002, PN0252_0003, PN0252_0004 are controls
and PN0252_0005, PN0252_0006, PN0252_0008 are cancers
I want to look at the difference of the proportion of cells in each cell type between cancer and control
But I am confused for getting proportion should I divide the number of cells on total cells in the object or total cells in each group
Basically I need a box plot like although here numbers are not correct

This looks like a ggplot with a facet added to it. What have you tried?

What I plotted here is what I have tried but I am not sure about the proportions

Your boxplot seems correct. What do you think is wrong about it?

I don't know proportion of T cells for example, I should divide the number of T cells in each patients on the total number of cells or total cells in controls and cancers separately

The answer to that depends on the question you're asking and what you wish to showcase. In any case, only you can answer it, we can't do that for you.

Ram RS as always your comments are extremely helpful . Thanks a million

@mahta-mira, are you able to write up your solution as an answer, then mark it as a selected answer? This will help the StackExchange bots work out what's going on.

At the moment I am not certain if my solution is correct :(


text:
Remove Redundant Sequences from FASTA file in Python
I'm attempting to remove redundant sequences from a fasta file (from NCBI). When I execute this code, it returns the number of spots, not the number of sequences. (Number of spots: 408,293, Number of sequences: 40,830)
Please let me know if you have run a code that can remove the duplicate ideas or have any ideas/suggestions.
Thank you SO much.
from Bio import SeqIO
import time

start = time.time()

seen = []
records = []

for record in SeqIO.parse("b4r2.fasta", "fasta"):
    if str(record.seq) not in seen:
        seen.append(str(record.seq))
        records.append(record)

#writing to a fasta file
SeqIO.write(records, "no_dupes_b4r2.fasta", "fasta")
end = time.time()

print(f"Run time is {(end- start)/60}")

Hi, your code was not formatting and it became a bit of a headache to format. Please invest some more effort into presenting your questions better.

Thank you so much, I'm just learning to use Stack Exchange! I appreciate it

Your code (as currently formatted) should generate the correct output, but it is not efficient. Make seen a set to improve these O(n) in calls to O(1) time complexity- this will speed it up a lot
seen = set()
records = []

for record in SeqIO.parse("b4r2.fasta", "fasta"):
    if record.seq not in seen:
        seen.add(record.seq)
        records.append(record)

Thank you for your input!


text:
Debugging manta 1.6.0 seg fault
Also posted on biostars
I'm trying to run manta on RHEL 8.  I built it using gcc 9.2.0 with --build-type=RelWithDebInfo
The following command throws a "segmentation fault, core dumped".  Running it under gdb, doing a "bt" when it crashes, consistently gets me a stack trace with 8733 frames, of which ~8700 are
searchRepeats (opt=..., index=aNumber, word=..., wordIndices=..., wordStack=..., repeatWords=...) at manta-1.6.0.release_src/src/c++/lib/assembly/IterativeAssembler.cpp:586

Yes, always that line.  aNumber ranges from 20 up to 33319.
manta/1.6.0/libexec/GenerateSVCandidates \
  --threads 10 --align-stats alignmentStats.xml\
   --graph-file svLocusGraph.bin --bin-index 0 --bin-count 1 \
  --max-edge-count 10 --min-candidate-sv-size 8 \
  --min-candidate-spanning-count 3 --min-scored-sv-size 50 \
  --ref GCA_000001405.15_GRCh38_no_alt_short_headers_nonACTG_to_N.fa \
  --candidate-output-file svHyGen/candidateSV.0000.vcf \
  --tumor-output-file svHyGen/tumorSV.0000.vcf \
  --chrom-depth chromDepth.txt --edge-runtime-log svHyGen/edgeRuntimeLog.0000.txt \
  --edge-stats-log results/stats/svCandidateGenerationStats.xml \
  --edge-stats-report results/stats/svCandidateGenerationStats.tsv \
  --tumor-align-file sample.bam

These are the top stack frames.  So I'm assuming the problem is running out of memory
#0  0x0000155554637b4e in _int_malloc () from /lib64/libc.so.6
No symbol table info available.
#1  0x0000155554639982 in malloc () from /lib64/libc.so.6
No symbol table info available.
#2  0x000000000070ddd5 in operator new (sz=81) at /research/bsi/tools/shared/gcc-9.2.0/libstdc++-v3/libsupc++/new_op.cc:50
        p = <optimized out>
#3  0x000000000078920d in std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >::_M_mutate (this=this@entry=0x15555419d1f0, __pos=__pos@entry=40, __len1=<optimized out>, __s=__s@entry=0x0, __len2=__len2@entry=1)
    at /research/bsi/tools/shared/build/x86_64-pc-linux-gnu/libstdc++-v3/include/bits/basic_string.tcc:310
        __how_much = 0
        __new_capacity = 80
        __r = <optimized out>

This is a loaded system.  How do I give manta more memory?

Has anyone seen this problem?

The bam file is only 23 gb

I ran it under  /usr/bin/time -v.  I got
Maximum resident set size (kbytes): 1084964

I've got a lot more than 1gb of ram available.  How do I get the manta tool to use it?

What is a loaded system? How are you checking that you have more than 1G RAM? Is it that it's running out of RAM (and, if so, do you have swap enabled?) or that a specific variable in the code is trying to store more data than the amount of memory allocated for that variable? Have you tried asking the manta devs?

https://github.com/Illumina/manta/blob/master/docs/userGuide/README.md#runtime-hardware-requirements mentions "<1GB per core" (emphasis mine). Please [edit] and tell us about your hardware.

@terdon Running on the Cluster I gave it 10 cores and 80GB.  Debugging is on a node that has 64 cores and a terabyte of memory.  It's doing a seg fault during a "new" call, so it's running out of "available to it" memory, but it's coming nowhere near to exhausting the memory available on the system

Just following up; did you work out the cause and/or a solution this problem @GregDougherty?

@jared_mamrot I built it with a much earlier version of gcc, and that fixed the problem

That's surprising! I tried to replicate the issue in a CentOS VM but failed, so I'm glad you figured out a solution :)

@jared_mamrot thank you.  I think it only happened to us once we got on RHEL 8

So, I "fixed" the problem by building with gcc 4.9.4, rather than gcc 9.2.0.
Why it's necessary to build with that much earlier, and completely out of date, version, I don't know.
But it worked


text:
SCTransform Warning: in theta.ml(y = y, mu = fit$fitted) : iteration limit reached
Running SCTransform on my seurat object produces the warning:
Warning in theta.ml(y = y, mu = fit$fitted) : iteration limit reached

What does this warning mean? Can I trust the output?
SCTransform runs through the data producing this warning until the progress bar reaches 100%. The resulting seurat object does have an SCT assay and I have successfully done downstream analysis with it. I am concerned that my SCTransform has not successfully done what it is supposed to do. Is there a way to cross-check?

The explanation can be found in link posted by ATpoint in his comment:

I get at least 10 posts on GitHub, GitMemory and Biostars just on the top page.

https://github.com/satijalab/seurat/issues/1378
https://github.com/satijalab/seurat/issues/2578
https://github.com/satijalab/seurat/issues/1378
https://www.biostars.org/p/450394/

The first author of the SCTransform paper, Christoph Hafemeister had posted this answer in a similar query on the satijalab GitHub page https://github.com/satijalab/seurat/issues/1378
"These warnings are related to the estimation of the theta parameter of the Negative Binomial regression model used to describe the molecule count distribution of each gene. When a gene has very few non-zero observations estimation of that parameter is not reliable and that's what the message indicates. This is not a problem, since the estimated theta parameters are not used directly for normalization but are regularized (sharing information across genes and thus lowering uncertainty)."


text:
Download all pubmed central article ids that have a keyword appearing in their titles/abstracts and also falling between two dates
I am using the following code to download all pubmed central article ids that have a keyword appearing in their titles/abstracts and also falling between two dates. The code seems to work when I select database = "pubmed" but when I switch to database = "pmc" I get empty list. For example, using the start and end dates as mentioned in the code I get 4 PubMed ids: ['36902449', '36880658', '36583446', '35239250'] and when I check manually, the first one among them (which is '36902449') is also present in PMC database, of course, the id is different (10003631) but same article. The code is expected to fetch that id. I'd appreciate it if someone could point out what I'm missing here.
import os
import datetime
from Bio import Entrez

# Set the email address to identify yourself to NCBI
Entrez.email = "example@gmail.com"

# Set the date range for the search
start_date = datetime.datetime(2023, 3, 1)
end_date = datetime.datetime(2023, 3, 22)

# Format the dates as strings in the correct format
start_date_str = start_date.strftime("%Y/%m/%d")
end_date_str = end_date.strftime("%Y/%m/%d")

# Build the search query with the date range
keyword = "dosage compensation"
# database = "pubmed"
database = "pmc"

# Use esearch to fetch the list of matching article IDs
if database == "pubmed":
    search_term = f'("{start_date_str}"[Date - Publication] : "{end_date_str}"[Date - Publication]) AND {keyword}[Title/Abstract]'
    handle = Entrez.esearch(db=database,
                            term=search_term,
                            retmax=100,
                            api_key = os.environ['ENTREZ_API_KEY'])
elif database == "pmc":
    search_term = f'("{start_date_str}"[Publication Date] : "{end_date_str}"[Publication Date]) AND {keyword}[Title/Abstract] AND open_access[filter]'
    handle = Entrez.esearch(db=database,
                        term=search_term,
                        retmax=1000,
                        api_key = os.environ['ENTREZ_API_KEY'])

record = Entrez.read(handle)
id_list = record['IdList']
handle.close()

# Print the number of matching articles
print(f"Found {len(id_list)} articles published between {start_date_str} and {end_date_str} in {database}")
# print(id_list)

There seems to be a problem with using the [Title/Abstract] field. What you have looks correct and should produce the expected result, since it should capture articles where the search term appears in either the title or in the abstract, or in both. But it does not. Strangely, you get the expected result using {keyword}[Abstract].
A workaround could be to OR the fields with:
f'("{keyword}"[Title] OR "{keyword}"[Abstract])'

Perhaps something like:
search_term = ' AND '.join([
    f'("{start_date_str}"[Publication Date] : "{end_date_str}"[Publication Date])',
    f'("{keyword}"[Title] OR "{keyword}"[Abstract])',
    f'(open_access[filter])',
])


text:
Separate multiple sequence into different file, one sequence per file
I have a file contain multiple sequence, and I want to separate them by "gene:" into different file.
example:
example.fa
>KQK21959
ATGGTGCTCGACCCCGTCGTCGGGGGCTCGGCCTCCCGCAAGGTCGAGCCCCTAACCCTA

>KQK14546
ATGCCTCCGCCGGACGCGGCCACCTACGCTCCCACGCTGCCAATGGCGACGTACCTCCCA

How to get results like this? One sequence per file, and after > is gene id, not cds id
KQK21959.fa
>KQK21959
ATGGTGCTCGACCCCGTCGTCGGGGGCTCGGCCTCCCGCAAGGTCGAGCCCCTAACCCTA

KQK14546.fa
>KQK14546
ATGCCTCCGCCGGACGCGGCCACCTACGCTCCCACGCTGCCAATGGCGACGTACCTCCCA

There's an app for that! The exonerate package (see here) includes a nice fastaexplode command that does exactly what you're looking for and will be much, much faster than anything we give you here since it is written in C and not a scripting language. Just install exonerate (you can do that with sudo apt install exonerate on Debian-based systems) and then:
fastaexplode file.fa

This will create these two files (file.fa is the original input file from your question):
$ ls
file.fa  KQK14546.fa  KQK21959.fa

A simple Biopython solution:
from Bio import SeqIO

for record in SeqIO.parse("example.fa", "fasta"):
    SeqIO.write(record, record.id + ".fa", "fasta")

Note, if you have multiple sequences with the same id, only the last sequence will be kept in the final output files

Thank you! It works excellent.

There are many ways to do it. Here is one way to do in Python.
import os

with open("input.fasta") as f:
    for line in f:
        if line.startswith('>'):
            record_name = line.lstrip(">")
            outfilename = f"{(record_name.rstrip())}.fasta"
            if not os.path.exists(outfilename):
                out = open(outfilename, 'w+')
            else:
                print('File already exists!')
        out.write(line)

That really depends on the tools you are using. If anything will do you can use the following python snippet.
If you know that each sequence is on a single line and your file isnt gigantic use this
EDIT
Updated so the file isn't read to memory at once to address @terdon point. Although, the above solutions are probably much better.
with open( "seqs.fasta" ) as f:
    out = ""

    for line in f:
        if not line.strip(): continue
        if line.startswith( ">" ):
            if not type( out ) == str: out.close()
            gene = line[1:].strip().split()[0]
            out = open( f"{gene}.fa", "w" )
            out.write( line )
            continue
        out.write( line.strip() )

This can also be done using awk. This is harder to read, but probably faster than python:
awk 'sub(/^>/,""){
        FILE=$0".fa";
    print ">"$0 >FILE;
        next
    }{print > FILE}' file.fa.

This will define a new output file FILE at each header line, whose filename is the header minus the > character, suffixed with .fa. And all lines in the file are appended to the current output file.
EDIT: Changed to a faster version found on: https://crashcourse.housegordon.org/split-fasta-files.html#using-awk-instead-of-csplitsed

This worked for me today:
# Load required packages
BiocManager::install("Biostrings")

library(Biostrings)

# Read the multi-genome FASTA file
multi_genome_fasta <- readDNAStringSet("*.fasta")

# Extract individual genomes and save as separate FASTA files
for (i in 1:length(multi_genome_fasta)) {
    header <- names(multi_genome_fasta)[i]
    sequence <- as.character(multi_genome_fasta[i])

    # Combine header and sequence
    combined_data <- paste(">", header, sequence)

    # Save the combined data as a separate FASTA file
    file_name <- paste0(header, ".fasta")
    writeLines(combined_data, file_name)

    cat("Saved", header, "as", file_name, "\n")
}

cat("Separation completed.\n")


text:
How to adjust cutoffs in Mutect2
I'm currently performing some initial analyses on targeted Illumina sequencing data generated with DNA from diseased tissue. My hypothesis is that one particular gene might contain somatic mutations in this disease. I'm currently running Mutect2 in tumor-only mode, using default settings.
As I have only found very few mutations (this result might well be plausible), I'd like to adjust the Mutect2 cutoffs to increase the sensitivity and see whether the result changes. For instance, I'd like to reduce VAF, depth of coverage and QUAL. Unfortunately, I can't work out which of the many arguments have to be modified and how: https://gatk.broadinstitute.org/hc/en-us/articles/360037593851-Mutect2
Can you help me out? Thanks!
gatk Mutect2 \
-R $fastaFile \
-I $inputFile \
--germline-resource $germlineResource \
-O $resultDir/$sampleName.Mutect.vcf.gz

Mutect2 uses a statistical model based on the available reads. Unlike with some other tools, there is not threshold based on depth or VAF or QUAL. The main statistic that does govern emission is the LOD (Log ODds) which is the log10 of the odds ("ratio") that the variation being inspected is a tumor relative to the probability that it is a technical artifact. The default threshold for this is 3.0 (corresponding to a odds ratio of 1:1000).
Perhaps you'd like to see lower probability results, in which case you should lower the threshold:
Using --tumor-lod-to-emit 2.0 should give you many more results, many of them might be wrong. The LOD isn't well calibrated, so I'd take it with a grain of salt....validate the results you get.


text:
Determining significance of a variable in a glm model
How would one determine the significance of a variable in a glm model?
If I, for example, have a dataframe like the one seen below, how would I determine if the origin of the sample has a significant effect on the value? (this is the number of enzymes capable of degrading the substrate of that matter)
Substrate    variable value origin
cellulose       M09    8    free
mannan          M12    2    free
glycogen        M65    2    free
chitin          M87    4    free
cellulose       M90    2    isolate
manan           M78    1    isolate
glycogen        M21    4    isolate
chitin          M21    1    isolate

So far I have tried:
mcomp = glm.nb(value ~ origin, data = my_data)

summary(mcomp)
Deviance Residuals:
    Min       1Q   Median       3Q      Max
-0.9625  -0.9047  -0.9047   0.1212   3.5232

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)   -0.01657    0.06571  -0.252  0.80097
originisolate -0.21911    0.08180  -2.679  0.00739 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.3418) family taken to be 1)

    Null deviance: 2053.5  on 2679  degrees of freedom
Residual deviance: 2046.3  on 2678  degrees of freedom
AIC: 6517.5

Number of Fisher Scoring iterations: 1

              Theta:  0.3418
          Std. Err.:  0.0186

 2 x log-likelihood:  -6511.4590

So free becomes the intercept and then isolate if significantly different from that. Does this mean Origin has a significant effect on the value?
Would the better approach be to do the following?:
mcomp = glm.nb(value ~ origin + Substrate, data = comb_data)
summary(aov(mcomp))
              Df Sum Sq Mean Sq F value Pr(>F)
origin         1     23   22.55   6.612 0.0102 *
Substrate     44   1445   32.84   9.631 <2e-16 ***
Residuals   2634   8981    3.41
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’

This shows me that origin and substrate have an effect on value if I understand correctly?

Second viewing of the question from what I can see -0.22 as a coefficient of origin is a strong negative association, so yeah it has a major impact. Its not how I would have done it, but that looks to be the result.

First viewing,
I'm going to throw my hat in here. We don't know what 'origin' is about, anyway just throw everything, i.e. each substrate and the origin into the same regression calculation. Check for a low-residual and preferably do a Q-Q plot, transform your data it this doesn't look good.
The key and the thing you are missing is your regression weights, without that I couldn't say very much. If the regression weight is near zero for 'origin' then it has zero impact. If the regression weight of 'origin' is positively greater than everything else ... I assume there are skewed distributions of 'substrates' between the 'origins'. If the regression weight of 'origin' is negative but still greater than all other regression weights then it is adversely affecting the 'value' you are seeking.
I don't know the experiment, the biological system or really the 'substrate' assays, so I can't comment any further.
The two issues I have are:

Doing an ANOVA on the output of a regression analysis doesn't make much sense to me. It is not something I would do, nor something in ML or GLM I've encountered.
Are you doing pairwise substrate/origin calculations? I presue not, but just in case this not how GLM works.

oh I was under the impression aov(summary(glm_model)) was a pretty common practice as it nicely summarises down to a variable level? Where do I find the regression weights?

When transforming the data do I have to transform both variables or can I transform just one of them until they fit the same distribution? and given that I am using a negative binomial glm does a uniform distribution between the variables matter as much?
I am unsure if there is a skewed distribution of variables between origins, that is what I am looking to find out really.

hi @Michael, very nice points here, yes the term anova is quite misleading sometimes. Normally people associate associate with the classical analysis of variance. In more "modern" uses, it can refer to the breakdown of deviance in a model, https://stats.stackexchange.com/questions/289226/what-are-the-differences-between-anovas-and-glms

It is like a log likelihood ratio test essentially, the "F" statistic instead of Chi-sq is used.

@StupidWolf So it is fine to aov(summary(glm_model)) to get a nice overview of the model? And by pairwise substrate/origin calculations does you mean doing the model on a subsatrate of the dataframe for each substrate @michael?

There is no better method, it's a matter of what you want to test or what is your question.
Using the anova() or aov(), test the terms collectively. For example, in your example with Substrate, the null hypothesis is that the coefficients are all zero, meaning cellulose =0, mannan =0 , ....
If the question is, "do the isolate samples have a higher value than origin samples?", then you can use your first model, where free is set as the reference and you test whether the effect of isolate is non-zero. Likewise you can do this for substrate and set of them as your reference. You can also do other pairwise comparisons using this model.
If the question is, "does origin have a significant effect on value, after controlling for substrate?", then you can use your second model.

Would I control for substrate by forming the model as such: value ~ origin + substrate ?

yes you do. you are giving every substrate a meaning  and on top of that estimating the effect of origin

And does the order in which I give the variable to the model, i.e. origin before substrate, indicate some sort of predefined importance to the model?

it doesn't define the importance, you can put it in any order. However if they are correlated, most of the variance will go to the first. I don't think this will apply in your case

see https://cran.r-project.org/doc/contrib/Faraway-PRA.pdf blocking design. A good ieda to read this

Thank you I will look at this :)


text:
Subsetting a Seurat object based on colnames
I'm trying to subset my seurat object based on colnames. I have gone ahead and labeled each cluster and now I want to subset all the colnames that are in Cancer_human for human_colnames and all the barcodes that are not in Cancer_human for mouse_colnames but get an error.
human_colnames = colnames(scData[,scData$cell_ann == "Cancer_human"])
mouse_colnames = colnames(scData[,!scData$cell_ann == "Cancer_human"])

Error in `[.Seurat`(scData, , scData$cell_ann == "Cancer_human") : Incorrect number of logical values provided to subset cells

I'm filtering on colnames not rownames. I want all the barcodes from the Cancer_human.

The code could only make sense if the data is a square, equal number of rows and columns. At the moment you are getting index from row comparison, then using that index to subset columns. Also, please provide a reproducible example data for testing, dput(myData).

@zx8754 i agree for the reproducible example (consider the pbmc data set that comes with Seurat) but not the square argument (it would be correct for matrices) but Seurat objects contain multiple things and is implemented such that for a Seurat object so, so[, i] and so[[j]][i] both give information about cell i.

If you want the column names to be preserved, you can subset the column names as a vector:
human_colnames = colnames(scData)[scData$cell_ann == "Cancer_human"]
mouse_colnames = colnames(scData)[!scData$cell_ann == "Cancer_human"]

Did you tried the subset function?
# assuming that scData@meta.data looks like this
#         orig.ident  ...  cell_ann
# ATTCAA  sample1          Cancer_human
# ATTTAT  sample1          Cancer_human
# ...
# AAATCA  sampleX          Cancer_mouse
# ATTTAA  sampleX          Cancer_mouse

Idents(scData) <- "cell_ann"
human_subset <- subset(scData, idents = "Cancer_human")
mouse_subset <- subset(scData, idents = "Cancer_mouse")
mouse_subset <- subset(scData, idents = "Cancer_human", invert = TRUE) # alternative way

P.S. It is not clear what you are trying to do. You should add to your post how your data looks like.


text:
Encoding different DNAs in pySBOL3?
I want to encode two different DNAs using pySBOL3. First dna is an engineered region usually composed by a transcriptional unit (TU) or a set of TUs. The second is vector which is a basically a plasmid that contains a set of the previously described dna. This is for differentiate between engineered on different plasmids.
Now for dna I used as type: sbol3.SBO_DNA and role: sbol3.SO_ENGINEERED_REGION and for vector I used as type: sbol3.SBO_DNA and role: SO:0000155 (plasmid).
I would like to have feedback on this choice. For example the plasmid could also be an engineered plasmid SO:0000637 or circular plasmid SO:0002211 or just circular SO:0000988.
Here is my actual implementation:
for dna_id in params.get('dna'):
    #this needs to be handlded by vector
    dnas = Dna.object().filter(dna_id)
    dna = sbol3.Component(dnas.name(), sbol3.SBO_DNA)
    dna.roles = sbol3.SO_ENGINEERED_REGION
    dna.name = dnas.name
    dna.description = dnas.description
    if dnas == None:
        print('No dnas specified')

for vector_id in params.get('vector'):
    vectors =Vectors.object().filter(vector_id)
    vector = sbol3.Component(vectors.name(), sbol3.SBO_DNA)
    vector.roles = SO:0000155
    vector.name =vectors.name
    vector.description = vectors.description()

Your choices match the same choices that I've made when designated types for inserts and vectors.
The more specific plasmid types are reasonable to use as well, but choosing just "plasmid" is certainly safe. That's why I have preferred to use that in my code because (1) my code might sometimes be applied to a different sort of plasmid, (2) I don't currently have a purpose for making the finer distinction, and (3) circular DNA would be indicated by SO.circular (SO:0000988) in any case.
Side note: I would suggest using the tyto package so that your code can use names like SO.plasmid in your code instead of opaque symbols like SO:0000155.


text:
dbSNP GRch37 Download
I am still new to this.
I am trying to download the GRch37 dbSNP data from the NCBI dbSNP Browser using the "Send to:" button. But whenever I select
"Send to:" > "File" > "Tab" > "Default order"
The coordinates of my positions will be GRch38 based. Does anyone knows how can I download the dbSNP data with GRch37 coordinates?
Thanks
Ken

I see no option to specify the reference version using the NCBI dbSNP browser, it is set to download the latest built information. I think this is not the standard tool to download dbSNP data across many sites, but rather to browse it and visualize it for specific variation.
If you are only interested in getting info for a smaller range of positions, I recommend you to use the UCSC browser, go to Table Browser, -> assembly GRch37, -> group variation, -> track all SNPs, -> table snp51 -> have a look at "describe table schema" to see what will you get, select your genomic range, and click on get output.
Alternatively and if you have a lot of space available in a server, simply go to the NCBI dbSNP FTP site, and download the latest GRch37 dataset (GCF_000001405.25) with all the entries. Make sure you read the README for further info about the files and data.


text:
How to convert vcf to tbed/tfam and again convert it back to tbed/tfam to vcf using plink
When converting vcf to tbed/tfam format, I am getting error like
unrecognized flag<'--transpose'>

Code I used:
plink2 --vcf sample_vcf_file.vcf --recode --transpose --out outputfile

Correct command is
plink2 --vcf sample_vcf_file.vcf --recode transpose --out outputfile

vcftools also can convert vcf into tped/tfam

Yep, from plink manual: --transpose Deprecated. Use "--recode transpose"..


text:
How to highlight the specific peptide sequences after performing multiple sequence alignment for the fasta file?
I have the peptide sequences and fasta files separately. I first aligned the fasta files using msa package. After that I'm trying to highlight the peptide sequences in the multiple sequence alignment output. I couldn't find a way to do that. Any suggestions on how to do this would be very helpful.
My R script so far:
# load required packages
library(seqinr)
library(msa)
library(dplyr)

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

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

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

Hi and welcome. Could you please define what you "specific peptide sequences"? Perhaps also include the title of the question in the question.

@M__ I have some short peptide sequences belonging to these fasta sequences, i want to highlight these peptides positions to know their locations in the msa output

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

I have used jalview myself. But I want to do it in a R program


text:
Error when using awk command to trim sequence files
I am attempting to edit my fastq sequence files with an awk command from a published article. The sequence files are all contained in one directory where I am running the awk command from the terminal in a Linux environment. There are 4 fastq files in the directory
01_thu159_IX_RNA_001.fastq
01_thu163_IX_RNA_001.fastq
01_thu164_IX_RNA_001.fastq
01_thu190_IX_RNA_001.fastq

I would like to trim all these files at once with the awk command so I wrote this for loop:
for filename in ./*.fastq*
do
awk -v s=10 -v e= 40 ' {if (NR%2 == 0) print substr ($0, s+1, length ($0)-s-e); else print $0;} ' ./${filename} >  ./${filename}_trimmed
done

However, I am getting the following error:
awk: fatal: cannot open file ` {if (NR%2 == 0) print substr ($0, s+1, length ($0)-s-e); else print $0;} ' for reading (No such file or directory)

I then thought maybe it was an error with the for loop so I tried the command as per the instructions in the article
awk -v s=10 -v e= 40 '{if (NR%2 == 0) print substr ($0, s+1, length ($0)-s-e); else print $0;}' ./01_thu159_IX_RNA_001.fastq > ./Trim_01_thu159_IX_RNA_001.fastq

But this also produced the same error:
awk: fatal: cannot open file `{if (NR%2 == 0) print substr ($0, s+1, length ($0)-s-e); else print $0;}' for reading (No such file or directory)

I have checked the syntax for the awk command and it should be correct. Could anyone please give me advice as to what I am doing wrong and how to fix this issue?

You have a whitespace here:
e= 40

awk thinks that 40 is your rule and the rest of the line is a filename.

Thank you very much really appreciate it, I didn't notice that. Will give it a try now again.


text:
Where can I get hg38 Copy number Variation (CNV) file?
I'm wondering how can I get standard hg38 Copy number Variation (CNV) file from UCSC website.
Although, I got a link that is not working: http://hgdownload.cse.ucsc.edu/goldenPath/hg38/encodeDCC/wgEncodeCnvGtapLatest.bed.gz

Can you share a genome browser URL that shows CN? I can't seem to find any CN track on the UCSC genome browser

What does "standard hg38 Copy number Variation" mean? What exactly are you expecting? Frequencies of known CNVs in the population? A list of known human CNVs? Rare CNVs that are actually present in the reference genome? Please [edit] your question and clarify.

If you want databases of known CNVs, you can have a look at:

Database of Genomic Variants (DGV)
DECIPHER
NIH GTEx
Clinvar CNV
NCBI dbVar

You can also use VarSome, a free genomic variant search engine developed by the company I work for, which (among other things) shows structural variants from all of the sources above in a single genome browser. Click on the "Structural Variants" card, and then:

Thanks for all the responses,
I found one database that caters standard CNVs of hg38 as well as other versions. Although, I've to confirm which dataset will going to work for me. Please let me know if anyone has better thoughts or suggestions.
http://dgv.tcag.ca/dgv/app/downloads?ref=GRCh37/hg19

I'd look at NCBI's resources first as recommended by @terdon. This will provide the basis for further comparisons, such as Varsome, NIH, TCAG etc ...