text:
Moving UMI from tag into read name
I have a bam file with Unique Molecular Identifiers (UMIs) for each read present on the RX tag ( such as RX:Z:TGAGAAGGG), as expected by picard and fgbio tools.
However, several other tools (such as UMI-tools) expect the UMIs to be present at the end of the read name, such as (NB501943:73:HK333BGXC:3:12402:20313:16437_TGAGAAGGG).
Does anyone know of a tool to do this?
I have an awk solution which I'll self-answer with, but I don't know whether that is robust.

@ATpoint, unless I'm incorrect, UMI-tools can't do this directly. umi_tools extract moves sequences from a position in the FASTA sequence onto the read name, but in this case the reads have been added to a tag via fgbio.

My current solution is (awk is not a tool I'm very familiar with, so this may not be good awk code, there is presumably a way to get it to print the header lines unmodified.):
samtools view -H temp.sam > temp2.sam
samtools view temp.sam |
awk '{OFS = "\t"} {for(i=1;i<=NF;i++)  if ($i ~ /^RX:Z:/) {$1=$1"_"$i; gsub("RX:Z:","",$1); print}}' >> temp2.sam

This seems to work when I test it, but I don't know if I'm missing anything.


text:
Where did the BIND website go?
I was trying to reach the BIND (Biomolecular Interaction Network Database) but the website is down.
Where can I find something similar? Has this platform moved somewhere or what?

Answer from @heathobrien converted from comment:
looks like that link has been dead for years:
https://www.biostars.org/p/16763/
Answer from @jason-aller converted from comment:
The biostars article mentions this:
http://baderlab.org/BINDTranslation

source: Stack Exchange

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

sort: votes


text:
Nvidia Parabrick fq2bam pipeline error - No such file or directory
I am trying to implement Nvidia-parabrick fq2bam pipeline to process my WES dataset. The reason for opting this pipeline is the huge number of samples (1004) which take many months to process for somatic mutation identification. THe GPU accelerated pipeline will be more helpful to get somatic mutations in a quick manner. After pulling parabrick docker image from here, the following command works fine and return the options :

docker run -it nvcr.io/hpc/parabricks:v2.5.0 -h

But whenever I give the following command to run fq2bam pipeline, I am getting an error of No such file or directory:

docker run -it nvcr.io/hpc/parabricks:v2.5.0 fq2bam --ref ucsc.hg19.fasta --in-fq sample_1.fastq.gz sample_2.fastq.gz --knownSites 1000G_phase1.indels.hg19.sites.vcf.gz --out-bam sample_dups_gpu.bam out-recal-file sample_recal_gpu.txt --tmp-dir /tmp

The error is as follows:
Traceback (most recent call last):
File "/parabricks/run_pipeline.py", line 7, in <module>
sys.exit(PB.pb_main())
File "PB.pyx", line 846, in PB.pb_main
File "/parabricks/pbargs.py", line 1287, in getArgs
return PBRun(sys.argv)
File "/parabricks/pbargs.py", line 392, in __init__
self.runArgs = getattr(self, args.command)(argList)
File "/parabricks/pbargs.py", line 447, in fq2bam
self.addToParser(bsmb_parser_sysgroup, sysOptionGenerator().allOptions)
File "/parabricks/pbargs.py", line 50, in __init__
PBOption(category="sysOption", name="--num-gpus", default=GetNumGPUs(), typeName=int, helpStr="Number of GPUs to use for a run")
File "/parabricks/pbutils.py", line 48, in GetNumGPUs
output = subprocess.check_output(["nvidia-smi", "-L"], universal_newlines=True)
File "/usr/lib/python2.7/subprocess.py", line 567, in check_output
process = Popen(stdout=PIPE, *popenargs, **kwargs)
File "/usr/lib/python2.7/subprocess.py", line 711, in __init__
errread, errwrite)
File "/usr/lib/python2.7/subprocess.py", line 1343, in _execute_child
raise child_exception
OSError: [Errno 2] No such file or directory

As I am using the docker image as it is,it is not clear to me which file is missing. Kindly suggest. Thanks.

Just a tip - if you edit your title to make it more informative (e.g. Nvidia Parabrick error - No such file or directory) it helps people answer the questoin and also helps other people who have the same error.

Thanks. I have approved the edit in question title to make it more informative.

My query was answered by NVIDIA Parabrick developer team. The installation instruction are available at https://forums.developer.nvidia.com/uploads/short-url/t5uXsZgdMXx0G2wDcYamFl9sO2x.pdf.

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


text:
How to filter out partially similar strings from two lists and make one list?
I have following two lists.
List1
m64041_200717_231916/100664277/0_42364 id=2440
m64041_200717_231916/100008029/0_42112 id=244006
m64041_200717_231916/100665802/0_46854 id=2440067
m64041_200717_231916/100008044/0_60161 id=24400
m64041_200717_231916/100928073/0_11804 id=2440152
m64041_200717_231916/10094492/0_44157 id=2440199
m64041_200717_231916/100008373/4102_50601 id=2440234

List2
m64041_200717_231916/100664277/0_42364
m64041_200717_231916/100665802/0_46854
m64041_200717_231916/100928073/0_11804
m64041_200717_231916/10094492/0_44157

These are sequence read ids from a genome file. List1 contains all the ids from the List2 but with some id numbers (id=2440). So the sequence ids are only similar to some extent. I want to compare these two files and I need to filter out List2 sequence ids from List one (with id numbers) as below.
Output.
m64041_200717_231916/100664277/0_42364 id=2440
m64041_200717_231916/100665802/0_46854 id=2440067
m64041_200717_231916/100928073/0_11804 id=2440152
m64041_200717_231916/10094492/0_44157 id=2440199

I will be grateful if someone can help me with this issue.

It sounds easy but I don't really understand exactly what you want to compare. Can U expand?

For another example imagine there are two lists. List1.
UnitedStates id=12,
UnitedKingdom id=15,
Spain id=20,
France id=25,
List2
UnitedStates,
Spain
I want to filter out UnitedStates id=12 and Spain id=20 using the List2. List1 and List2 are two files.
Output should be,
UnitedStates id=12,
Spain id=20

ah, so only keep whats in list 2?

Yeah. list2 will be the query.

You really need to learn some python :)
with open('/path/list2.txt', 'r') as fin:
    to_keep = set([])
    for line in fin:
        to_keep.add(line.strip())

with open('/path/list1.txt', 'r') as fin:
    with open('/path/list3.txt', 'w') as fout:
        for line in fin:
            if line.strip().split(' ')[0] in to_keep:
                fout.write(line)

Thank you very much. It worked. Yes, I need to learn python.

And another question. Do you know a way to filter out/remove those sequences from the original fasta file using the id list created from the above code? I have a piece of code I used earlier to remove sequences using a list. But since these ids contains a space, It's not working now. I tried to change it, but it didn't work. It'll be an easy thing for you,but I couldn't figure it out. Can you please help?

Sure. Post a new question with an example plz

Why would you hardcode two file names instead of just reading argv for the input file and printing to stdout so you can then redirect?

He asked for files

This can be done in python easily. If you have different text files then you can load them into a python list and follow from step 2
# Step 1
list1 = ["m64041_200717_231916/100664277/0_42364 id=2440",
"m64041_200717_231916/100008029/0_42112 id=244006",
"m64041_200717_231916/100665802/0_46854 id=2440067",
"m64041_200717_231916/100008044/0_60161 id=24400",
"m64041_200717_231916/100928073/0_11804 id=2440152",
"m64041_200717_231916/10094492/0_44157 id=2440199",
"m64041_200717_231916/100008373/4102_50601 id=2440234"]

list2 = ["m64041_200717_231916/100664277/0_42364",
"m64041_200717_231916/100665802/0_46854",
"m64041_200717_231916/100928073/0_11804",
"m64041_200717_231916/10094492/0_44157"]

# Step 2
split_list1 = [ele.split(" ")[0] for ele in list1 ]
mapdict = {split_list1[i]:list1[i] for i in range(len(list1))}
common = [value for value in split_list1 if value in list2]
required = [mapdict[x] for x in common]
print(required)


text:
How to get rows with similar values in two different columns using command line?
I have following example data from a large file.
cont100 1128 1125
cont1005 3642 3642
cont1006 103 19
cont1037 3146 3146
cont104 895 890
cont1056 934 934
cont1059 1750 1750
cont1072 2577 2577
cont1078 43 42
cont1081 727 2

I need to get only the rows with similar values in column 2 and column 3. The output should be like this.
cont1005 3642 3642
cont1037 3146 3146
cont1056 934 934
cont1059 1750 1750
cont1072 2577 2577

Also I need to get the inverse of that. That means I need to get the rows with values that are not similar in the columns. The output should be like this.
cont100 1128 1125
cont1006 103 19
cont104 895 890
cont1078 43 42
cont1081 727 2

Do you mean identical values in columns 2 and 3?

I think it is easy to do with awk, e.g. see here.
For your example, with your data in file:
%  awk '$2 != $3 {print $0}' file
# outputs
cont100 1128 1125
cont1006 103 19
cont104 895 890
cont1078 43 42
cont1081 727 2

%  awk '$2 == $3 {print $0}' file
# outputs:
cont1005 3642 3642
cont1037 3146 3146
cont1056 934 934
cont1059 1750 1750
cont1072 2577 2577

Note that where you say "similar" I am interpreting that to mean "exactly the same", because that is the case in your example.

It worked and I got the result that I wanted. Thank you very much.

The Perl solution is similar to the awk solution from MaximilianPress:
Print lines with identical values in columns 2 and 3:
perl -lane 'print if $F[1] == $F[2];' in.txt

# Prints:
cont1005 3642 3642
cont1037 3146 3146
cont1056 934 934
cont1059 1750 1750
cont1072 2577 2577

Print lines with different values in columns 2 and 3:
perl -lane 'print if $F[1] != $F[2];' in.txt

# Prints:
cont100 1128 1125
cont1006 103 19
cont104 895 890
cont1078 43 42
cont1081 727 2

The Perl one-liner uses these command line flags:
-e : Tells Perl to look for code in-line, instead of in a file.
-n : Loop over the input one line at a time, assigning it to $_ by default.
-l : Strip the input line separator ("\n" on *NIX by default) before executing the code in-line, and append it when printing.
-a : Split $_ into array @F on whitespace or on the regex specified in -F option. Arrays in Perl are 0-indexed.
SEE ALSO:
perldoc perlrun: how to execute the Perl interpreter: command line switches


text:
Association analysis for methylation beta value and phenotypes in R
I am trying to learn how to do association analysis of methylation array data with numerical trait like body weight, eGFR.
As I have seen some of the publications, methylation data association with eGFR has been showed through regression analysis but I am not sure how to run regression for so many CpG probes in R.
Can you have look on this below publication and could you help with R code for regression analysis?
https://www.nature.com/articles/s41467-017-01297-7.pdf
I do have methylation beta value from RnBeads but I am not sure how to read methylation data (rows/column-CpG/SampleID) and then how to read phenotype data file and finally how to make single sheet to run regression analysis?

Have a look at this method https://academic.oup.com/bioinformatics/article/34/15/2657/4939328

Hi, I found a way to do methylation data association analysis with eGFR as shown here (https://support.bioconductor.org/p/93868/), but not sure how make input files; like in methylation data file; do I take CpG probe as rows or column? and I am also not sure how to mention covarites in modal.matrix for eGFR association? Many thanks

Answer from @user4939328, converted from comment. Please improve this answer if you can:
Have a look at this method:
https://academic.oup.com/bioinformatics/article/34/15/2657/4939328
[Wondering why this was converted from a comment? See the StackExchange rules about comments here]


text:
How to remove sequences from a fasta file using a sequence ID list which contains a space within the id?
I have a fasta file that contains sequence reads and sequence id file that needed to be removed from the fasta file. I have done this earlier, but since id contains a space my piece of code is not working and I don't know how to change it since I'm still learning. fasta file looks like this.
>m64041_200717_231916/74/12723_24868 id=30
CCGCACCTCCTCAATCTGCAGCAGTTGAGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCA
>m64041_200717_231916/77/1941_50622 id=3115
TGAGCCGCACCTCCTCAATCTGCAGCAGTTGAGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCATGAG
>m64041_200717_231916/105/20691_65844 id=488
AGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCAAGACTTTAT
>m64041_200717_231916/108/17414_67048 id=4956
TGAGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCATGAGGCAGACTTTATCATAGTCACTCACACTTG
>m64041_200717_231916/162/0_6615 id=857
CAGACTTTATCATAGTCACTCACACTTGTCCATGAGGCAGACTTTATCATAGTCACTCACACTTG

IDList file is like this.
m64041_200717_231916/74/12723_24868 id=30
m64041_200717_231916/108/17414_67048 id=4956
m64041_200717_231916/105/20691_65844 id=488

I need to remove these sequences from that fasta file. I have a code I used earlier. But it is not working since id contains a space (I think). Previous code is like this.
# use as follows
# ./filter_seq.bash <fastafile> <list_of_ids_you_want_to_remove>
# ./filter_seq.bash all.fasta bacterial_IDlist

grep ">" $1 | sed 's/>//g' | sort | uniq > all_del
sort $2 | uniq > query_del
comm -23 all_del query_del | sort | uniq > filtered_del
perl -ne 'if(/^>(\S+)/){$c=$i{$1}}$c?print:chomp;$i{$_}=1 if @ARGV' filtered_del $1 > filtered.fasta
rm *_del

I would be grateful if someone can help me.

I would highly recommend getting rid of the spaces first using something like sed 's, ,_,g' -I in.fa. Spaces in any kind of genetic data file are bad news and will cause downstream errors. It will make further analysis much easier.

@user438383 why do you say that? Spaces in fasta headers are incredibly common, allowed by the standard and should not cause any issues at all. I can't think of any tool that would have trouble parsing a fasta header like that and would even say that if a tool cannot do this, it is not fit for purpose. As an example, have a look at pretty much any fasta sequence from refseq. Here's human tp53: https://www.ncbi.nlm.nih.gov/nuccore/NG_017013.2?report=fasta.

Install (pip) biopython if u don't have it.
from Bio import SeqIO
with open('/path/ids.txt', 'r') as fin:
    to_keep = set([])
    for line in fin:
        to_keep.add(line.strip().split(' ')[0])
with open('/path/output.fa', 'w') as fout:
    for record in SeqIO.parse('/path/sequences.fa','fasta'):
        if record.id in to_keep:
            SeqIO.write(record, fout, 'fasta')

I tend to do this using the FastaToTbl and TblToFasta scripts I have posted before. FastaToTbl will convert fastq sequences to single-line, tab-separated recors which makes them far easier to parse. Using these two scripts, you can do:
$ FastaToTbl file.fa | grep -vwf idList | TblToFasta
>m64041_200717_231916/77/1941_50622 id=3115
TGAGCCGCACCTCCTCAATCTGCAGCAGTTGAGGCCACTACCCTTCTGCTCAATGGTTCC
TGCAGACTTTATCATAGTCACTCACACTTGTCCATGAG
>m64041_200717_231916/162/0_6615 id=857
CAGACTTTATCATAGTCACTCACACTTGTCCATGAGGCAGACTTTATCATAGTCACTCAC
ACTTG

The -v option will tell grep to return non-matching lines, and the -f option tells it to read the patterns to search for from a file. In this case, we give it the file with the IDs to skip (idList). The -w ensures we only find "whole word" matches, so that something like m64041_200717_231916/77/1941_50622 id=31151111 will not count as a match for m64041_200717_231916/77/1941_50622 id=31151111.

Another approach would be to just use awk directly:
$ awk '(NR==FNR){s[">"$0]++} (/^>/){ s[$0] ? a=0 : a=1}a' idList file.fa
>m64041_200717_231916/77/1941_50622 id=3115
TGAGCCGCACCTCCTCAATCTGCAGCAGTTGAGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCATGAG
>m64041_200717_231916/162/0_6615 id=857
CAGACTTTATCATAGTCACTCACACTTGTCCATGAGGCAGACTTTATCATAGTCACTCACACTTG

Here is the same awk script slightly expanded to make it easier to understand:
 awk '{
        if (NR==FNR){idsToSkip[">"$0]++}
        else if (/^>/){
            if($0 in idsToSkip){
                wantThisId = 0
            }
            else{
                wantThisId = 1
            }
        }
        if(wantThisId == 1){
            print
        }
      }' idList file.fa


text:
What reasons are there to choose Illumina if PacBio provides longer and better reads?
PacBio provides longer read length than Illumina's short-length reads. Longer reads offer better opportunity for genome assembly, structural variant calling. It is not worse than short reads for calling SNP/indels, quantifying transcripts. Sounds like PacBio can do whatever Illumina platform can offer.
Would that be any reason for going with Illumina if PacBio's long reads can do everything and more?

There are so many reasons why one might want to prefer Illumina over PacBio (also note that it's a false dichotomy, at least Oxford Nanopore is a competitive sequencing platform):

The first (IMHO and the most common reason) is still the cost of both sequencing and the instruments. Illumina can sequence a Gbp of data for \$7 - \$93. PacBio sequencing is according the same webpage \$115 per Gbp, however at our sequencing center it's ~$200. Though ONT might have already a cheaper solution. Edit, I just found a google sheat with prices that seems to be frequently updated, seems that the ratio still holds Illumina short reads ~10x cheaper than PacBio.
RNA-seq (i.e. analysis of a gene expression) is not possible with PacBio due to preferential sequencing of smaller fragments; shorter genes would always be shown to be more expressed. To be clear, it's possible to sequence RNA with PacBio (the keyword is iso-seq), but the analysis of gene expression is problematic.
It's way easier to extract fragmented DNA (concerns small non-model organisms; although recently a single mosquito was sequenced, so we can expect a further improvement)
other sequencing techniques as RAD-seq that allow genotyping with very little effort and cost, I have never seen anybody even considering using long reads for such genotyping
Genome profiling (assembly-less genome evaluation) based on kmer spectra analysis is not possible with PacBio data due to higher error rates. Conflict of interest: I am a developer of one of the tools for genome profiling (smudgeplot)
DNA in Formalin-Fixed Paraffin-Embedded samples is fragmented anyway (~300bp), therefore, there is no point in sequencing them with more expensive long-read technology (contributed by @mRoten)
I bet there will be a plenty other applications I am not aware of

Please add that FFPE sample DNA is fragmented (~300 nt), so Illumina short read sequencing is still has the lowest cost per base than PacBio. Even ConcatSeq (https://www.nature.com/articles/s41598-017-05503-w) only increase the read count by a factor of 5. Illumina is still the winner.

I turned the answer to community wiki so it's easier for people to add/updates reasons for Illumina. As this is not complete not a stable-in-time answer.

Many analyses performed on Illumina machines these days require large numbers of reads. For example, most analyses in ChIP-seq, RNA-seq, ATAC-seq etc, need 10s or even 100s of millions of reads for the statistics to work out properly.
But this isn't limited to just sequencing as an assay experiments. High depths are important for things like somatic variant calling as well. Or simply sequencing of medical gene pannels, where you might only have 20 amplicons of a 500bp each, but need to do it on 100s of patients as cheaply as possible.

Three reasons for Illumina:

Much better for a large number of samples (easily handle 96 samples).
SNP calling is much better - much greater depth
Hardware costs, an Illumina MiSeq machine is cheap

PacBio SNP calling has been improved towards the standard of Illumina by DNA modification to the inputs to be sequenced but it depends on the application.

As others have mentioned, cost is probably one of the main issue.
Other than that, PacBio sequencing also requires more input DNA material (~3-10 ug for human) while Illumina sequencing works with just 100ng. This makes PacBio sequencing more challenging when samples are rare. I'm sure people at PacBio are trying to resolve some of these concerns already though.


text:
Different results in differential expression analysis of microarray data
I am performing differential gene expression analysis to microarray data for type 2 diabetes donors and nondiabetic donors. When I run the code I get some different results in each time (about 50 or so different genes in each run). Is this normal? I am submitting a project and it could be very problematic if my results change in each run, since I will perform downstream analysis for the DE genes and use web tools such as DAVID for further analysis.
Is there a way to keep the results from changing in each run?
g <- getGEO("GSE20966",GSEMatrix = TRUE)
g<- g[[1]]
exprs(g)
exprs(g) <- log2(exprs(g))

#Diffrential gene exp

library(limma)
design <- model.matrix(~0+sampleInfo$characteristics_ch1.1)

colnames(design) <- c("NonDiabetic","T2D")

summary(exprs(g))

exp <- na.omit(exprs(g))
cutoff <- median(exp)

## TRUE or FALSE for whether each gene is "expressed" in each sample
is_expressed <- exprs(g) > cutoff

## Identify genes expressed in more than 2 samples

keep <- rowSums(is_expressed) > 2

## check how many genes are removed / retained.
table(keep)
keep<- na.omit(keep)

## subset to just those expressed genes
g<- g[keep,]
fit <- lmFit(exp, design)

head(fit$coefficients)
contrasts <- makeContrasts(T2D - NonDiabetic, levels=design)
fit2 <- contrasts.fit(fit, contrasts)
fit2 <- eBayes(fit2)
topTable(fit2)
topTable(fit2, coef=1)
decideTests(fit2)
table(decideTests(fit2))
aw <- arrayWeights(exp,design)
aw
fit <- lmFit(exp, design,
            weights = aw)
contrasts <- makeContrasts(T2D - NonDiabetic, levels=design)
fit2 <- contrasts.fit(fit, contrasts)
fit2 <- eBayes(fit2)
g <- na.omit(g)
anno <- fData(g)
anno
anno <- select(anno,`Gene Symbol`,ENTREZ_GENE_ID)
fit2$genes <- anno
topTable(fit2)
full_results <- topTable(fit2, number=Inf)
full_results <- tibble::rownames_to_column(full_results,"ID")
library(ggplot2)
dev.off()
ggplot(full_results,aes(x = logFC, y=B)) + geom_point()
p_cutoff <- 0.05
fc_cutoff <- 1

full_results %>%
  mutate(Significant = adj.P.Val < p_cutoff, abs(logFC) > fc_cutoff )
 %>%
  mutate(Rank = 1:n(), Label = ifelse(Rank < topN, anno$`Gene
 Symbol`,"")) %>%
  ggplot(aes(x = logFC, y = B, col=Significant,label=Label)) +
geom_point() + geom_text_repel(col="black")

 p_cutoff <- 0.05
 fc_cutoff <- 1

filter(full_results, adj.P.Val < 0.05, abs(logFC) > 1)
 library(pheatmap)
topN <- 20

top<-150
gene_names <- mutate(full_results, Rank = 1:n()) %>%
  filter(Rank < top) %>%
  pull(Gene.Symbol)

DE results do not change from one run to another. The code given in your question will give identical results each time you run it.


text:
Why does gene count increase drastically after scaffolding with Hi-C data?
I have a conceptual question.
I have a diploid, outcrossing plant genome assembly of ~1.1 Gb size. The original assembly is generated from PacBio reads. After genome annotation with the Mercator software, the original assembly returned a gene count of 59,615. Then, the diploid genome assembly was phased with Hifiasm+Hi-C data. There, after running the annotation on two haplotype assemblies, the gene count was 75,936 for haplotype 1 and 67,233 for haplotype 2.
Then, we sent the pacbio assembly together with the haplotypes to Dovetail genomics for scaffolding. After the scaffolding with HiRise and annotation with Mercator, I got the following gene count for the three assemblies:
Scaffolded_assembly- 75,874
Scaffolded_haplotype_1- 67,243
Scaffolded_haplotype_2- 86,465

Can anyone please shed some light on why there is such an improvement in gene count? Also, the gene count for haplotypes seems to be the opposite after scaffolding. Do you have any ideas? I do not know how many of those are high confidence or low confidence gene or splice variants.
The BUSCO score was similar in all the cases.
It would be nice if you can give some feedback or any paper explaining such phenomena. Thanks.


text:
Number of collinear genes are zero while running MCScanX
I am trying to do a synteny analysis using MCScanX. As per requirement I have done all by all blast between two sets of protein from two closely related species. And also prepared the gff file. The file formats are as follows:
cyp_dan.blast
NP_001347641.1  XP_042591522.1  78.31   166     32      1       1       166     1       162     5e-91   267
NP_001347641.1  XP_042577882.1  60.24   166     62      2       1       166     1       162     4e-63   196
.
.
XP_704272.2     XP_018953536.1  31.03   87      59      1       275     360     90      176     8e-07   50.4
XP_704272.2     NP_001005958.2  33.33   63      38      1       273     331     38      100     1e-06   48.9
XP_704272.2     XP_018969808.1  31.33   83      56      1       275     356     91      173     2e-06   49.3```

cyp_dan.gff
NC_056572.1     gene-LOC1221465.1       2966    7093
NC_056572.1     id-LOC122146575.1       2966    3166
NC_056572.1     id-LOC122146575.1       3246    3357
.
.
NC_002333.2     gene-trnP;Dbxre.1       16527   16596
NC_002333.2     rna-trnP;Parent.1       16527   16596
NC_002333.2     exon-trnP-1;Par.1       16527   16596

cmdline: /home/software/MCScanX/MCScanX cyp_dan
Output:
############### Parameters ###############
# MATCH_SCORE: 50
# MATCH_SIZE: 5
# GAP_PENALTY: -1
# OVERLAP_WINDOW: 5
# E_VALUE: 1e-05
# MAX GAPS: 25
############### Statistics ###############
# Number of collinear genes: 0, Percentage: 0.00
# Number of all genes: 108556
##########################################

Here, both of the species are from same family since there is no chance to have "zero collinear gene".

I'm not familiar with the MCScanX tool, but I see two potential issues that are worth investigating.

The cyp_dan.gff file does not look like GFF or any of its common variants (GFF3, GTF, etc.). Again, I'm unfamiliar with MCScanX, so I don't know precisely how its wants the input data formatted. But if it wants a GFF file, I would double-check to make sure that what you have matches what the tools expects/requires.

The gene/protein IDs in the cyp_dan.gff file do not seem to match the query or subject IDs in the cyp_dan.blast file. No tool is going to be able to find collinear adjacent genes if the IDs in the BLAST results don't match the IDs in the gene annotations.

Hi,
Thanks for your response.

Original gff format is not allowed to use. It requires modification with having chr name, geneID, start and end position only.
I took out the necessary info using following cmd:
awk -F "\t" '{a=substr($9,4,15)}$0~/gene/{print $1 "\t" a".1" "\t" $4 "\t" $5}' A.gff > A_new.gff

I have extracted the info from the NCBI annotated file and didnt make any change in ID info by my self. And after extracting this is how it looks. Can you please give me any suggestion to find exact match between blast and gff file?


text:
Why isn’t ace2 included in SARS‑CoV‑2 PPI map?
I'm doing my master's thesis on drug repositioning for SARS‑CoV‑2 using AI more specifically network medicine,
I found only one human protein-protein interactions (PPIs) map from Gordon et al.
The manuscript doesn't include the ace2 protein which the SARS‑CoV‑2 spike protein directly binds with. AVE2 is a host cell surface receptor facilitating virus entry and thereafter replication.
Thus, I'm looking for PPI map that includes ACE2 protein but i can't find any
Secondly, is there a reason why it's not been included ?
PS: I have very basic knowledge of this domain

In fairness to the authors they do included it, viz. "captopril  modulates  ACE2  indirectly  via  its  direct  interaction  with  Angiotensin Converting  Enzyme,  ACE". So they hold results on ACE and therein have a claim for ACE2


text:
Where to find publicly available Sanger chromatography data?
I am looking for any public available databases where I can download Sanger chromatography data ideally in Ab1 or SCF file format. I need extremely large amounts of this, since I want to use it for machine learning.

@ATpoint , can you link me here any ML based solution for basecalling if you know any ? I found only 2-3 papers in the literature

Why do you want to reinvent the wheel on this? This is an old technology that doesn't require any new tools, not least of which because its use is decreasing every year.

@DevonRyan I can fully unterstand what do mean but Sanger DNA sequencing is not going away any time soon, also we are looking for inhouse Basecalling method/software which will be independent of any licencing

I must agree with the others. Furthermore, I'd say writing a script that can deal with frameshifted sequences in the PCR pool (the major issue with basecalling) is easier, but that is application specific as in general Sanger is used for validation and the data has to be perfect. E.g. I wrote a tool to call secondary peaks to test the diversity of error prone PCR pools, which is a super specific application.

There are some interesting trends in a Sanger trace, i.e. peak height distribution, which tell you about how well the experiment was (concentrations, gunk etc.), but that again is not too interesting to be heavily used —unless say implemented in Benchling or other major tool, i.e. the tricky part of a new tool is getting to be adopted.
Furthermore, about licensing the Ab1 files are what is spat out by the machine and that is where the licensing/compression comes into play, not in basecalling as one can very easily do the basecalling from the raw data within them ignoring the calls.

NCBI Trace Archive
The NCBI Trace Archive is a permanent repository of DNA sequence chromatograms (traces), base calls, and quality estimates for single-pass reads from various large-scale sequencing projects. Contaning more than 2 billion traces.


text:
Can an alternate allele be more common than a reference allele?
I'm trying to analyze data from a genotype-phenotype association study. The genotype data is stored in csv files that record for each variant site and each subject, the alleles that the subject has at that site. A quick online search reveals that 0 means the reference allele, 1 means the most common alternate allele, 2 means the 2nd most common alternate allele, etc.
But I noticed that at some variant sites, most (>90%) subjects have 1/1, i.e. two copies of the alternate allele. I asked around and some suggested I just swap the roles of reference and alternate alleles whenever I find that the alternate allele frequency > 0.5. My questions are:

Can there be a good reason to classify a less common allele as the reference allele? I thought that the reference allele was defined as the most common one in the population.
Is swapping the roles of reference and alternate alleles whenever I find that the alternate allele frequency > 0.5 a standard practice?

Thanks

You are mixing up a bit the terminology here. One thing is a reference allele and another is a major allele. The first refers to the genotype present at a given position in the "reference" genome (e.g. if you work with human data this can be the reference version hg19 or hg38, etc). The second refers to the most common allele in a population. This definition can lead to problems since what is the most common allele can differ among different populations. In this respect, it is perfectly fine to have a reference allele with a lower frequency than another alternate allele.

No for the reasons above. If you mean whether swapping major to minor is standard practice whenever you find out the major allele is actually minor as per the definition (less frequent), then I guess you can but this will lead to inconsistencies whenever you try to merge your data with other data sets, or other people merge their data with your data. Since this for sure happens and issues about strandness are well known, one has to make sure that you make explicit how you treat your data once you release it, and also one has to make sure that major/minor are defined the same way when merging data, and do the fix whenever needed.

Yes, depending on how the reference allele was determined. If, for example, you want to refer to hg38 as the reference genome, then it's important to use the hg38 allele as reference.

Without stating the reference allele convention? Not that I'm aware of. However many people refer to minor/major alleles rather than non-reference alleles, in which case it is assumed that the major allele is the most common allele.

As long as you state what your reference genome and/or population is, then it's fine to choose your own convention about reference alleles. But please be consistent and explicit about what you've done. One of my biggest struggles in GWAS has been working out the correct strand orientation for stated variants (important when merging datasets), and it doesn't help for complementary variants (i.e. A/T, C/G) when the ref/alt alleles are arbitrarily chosen.


text:
Create GFF3 feature exporter
Problem Statement

The GFF3 format is a commonly-used one in bioinformatics for
representing sequence  annotation. You can find the specification
here:
http://www.sequenceontology.org/gff3.shtml
The genome and annotation for Saccharomyces cerevisiae S288C is on the
class server here:
/home/jorvis1/Saccharomyces_cerevisiae_S288C.annotation.gff

Note that this same file has both the annotation feature table and the
FASTA sequence for the  molecules referenced. (See the '##FASTA'
directive in the specification.) Within the feature table another
column of note is the 9 th , where we can store any key=value  pairs
relevant to that row's feature such as ID, Ontology_term or Note. Your
task is to write a GFF3 feature exporter. A user should be able to run
your script like this:
$ export_gff3_feature.py --source_gff=/path/to/some.gff3 --type=gene --attribute=ID --value=YAR003W

There are 4 arguments here that correspond to values in the GFF3
columns. In this case, your  script should read the path to a GFF3
file, find any gene (column 3) which has an ID=YAR003W (column 9).
When it finds this, it should use the coordinates for that feature
(columns 4, 5 and 7)  and the FASTA sequence at the end of the
document to return its FASTA sequence. Your script should work
regardless of the parameter values passed, warning the user if no
features were found that matched their query. (It should also check
and warn if more than one  feature matches the query.)
The output should just be printed on STDOUT (no writing to a file is
necessary.) It should have a header which matches their query, like
this:

gene:ID:YAR003W .... sequence here ...

As an extra challenge, you can format the sequence portion of the
FASTA output as 60-characters per line, which follows the standard.

So far I have this below. I'm not even sure if this is correct to start with:
#!/usr/bin/env python3

for line in open("Saccharomyces_cerevisiae_S288C.annotation.gff"):
   line1 = line.rstrip()
   if line.startswith("#"):
       continue
   column = line1.split("\t")
   if len(column) != 9:
       continue
   id = column[8]
   type = column[2]

Can anyone help me with this? I'm totally lost. Thank you!

Was there any restriction about using already developed packages for this purpose? If not : https://biopython.org/wiki/GFF_Parsing , you can wrap your script around with this plus argparse to work from the console: https://docs.python.org/3/library/argparse.html

@Avamys please add this as a separate answer. Comments should only be used for transient information that can help improve or supplement questions or answers. See here for more information.

Your answer looks like a good start. You'll probably want to use sys.argv, argparse, or another argument parsing tool to read in options from the command line, then compare them to the fields (e.g. the "type" argument gets compared with the type column).
Given that attributes can be arranged in an arbitrary order, you'll probably need to read all of them from each line and store in a dict in order to determine whether the target is present. A regular expression search is another way to do that.
There are a few libraries available (e.g. from biopython) for reading and subsetting FASTA records.


text:
How to determine a correct promoter sequence? (Bisulfite sequencing)
I'd like to ask you for help with an university assignment. First let me post the entire question:

Design primers for bisulfite sequencing of the promoter of the ACE2
gene in microbat  Myotis lucifugus (based on the annotation of the
Broad Institute Myoluc2.0/myoLuc2). The gene is localized on the
contig GL429816. As a transcription start, use the start of the ACE2
transcript variant XM_023753670.1; once you find ACE2 in the UCSC
genome browser .You will find it in NCBI RefSeq track

I don't want anyone to solve it for me but I do have a problem with one aspect. I am not clear on how to determine what is a valid promoter for this gene in microbat. I do know it is the upstream sequence however I am not clear on how to tell the number of upstream bases that are the actual promoter. I did find some sources for the actual promoter sequence determined by much smarter people than myself however those are for ACE in humans so I am afraid that it is not applicable in my case.
I'd like someone to point me in the right direction where I could find my answer would be enough.

Thank you I decided to use 900bp, just because that is the legnth the teacher used in of his lectures. I would like to select you comment as an answer, but it's not possible it seems.


text:
What tools can I use to look up my alleles, genotypes, or phenotypes from my sequenced DNA (WGS)?
I bought a "Whole Genome Sequence" kit, which gave me the following VCF files (this service did not include any analysis, just the raw data):
cnv.vcf
indel.vcf
snp.vcf
sv.vcf
My understanding is that these files show the differences of my DNA relative to a reference human genome (hg19 in my case)
What tools and resources can I use to find out which known alleles/variants of genes I have?
An example: Blood type
According to Wikipedia, there are 6 common alleles of the ABO gene, and which alleles I have determine my blood type.
Looking at my VCF results, I have these SNPs and INDELs within that gene:
<removed>

Given this, how would I figure out my blood type?
I understand that my data may not be 100% accurate, and that many phenotypes are the result of more than one gene, environmental factors, or are not well studied at all, but this should work for some simple things ie blood type, COMT, BHLHE41, CYP1A2, etc

For the blood locus. Try and find online the genomic position (i.e. the number like 133273540 in the POS field of the .vcf) of the ABO blood locus and look it up on your dataset. Remember to make sure you find it on the gr37/hg19 build.

Wait, did you got yourself yourself sequenced?? I thought "my DNA" means "my DNA sample" not "my own DNA"!

And do I get right, you know your variants and you would like to know something about them, right?

Please don't just share your own genetic information online, especially if you are not fully understanding what you are looking at. It is not something to be posted everywhere, although some loci might be relatively harmless you shouldn't take this lightly.

Please ping me if you decide you want this censored. You should be very aware of the long-term consequences for you and your current and future relatives of sharing your genetic information on the internet in a post.

@DevonRyan thanks, I've removed it. I didn't think a small snippet would be a problem

@WouterDeCoster also thanks for the warning

If this is your own DNA and you are trying to analyse yourself you should talk to a genetic counselor. I personally think it's cool, but I deeply understand the potential consequences.
What you want is a variant annotation tool. There are many, but I normally write my own. Which is simply a case of pulling down a database and using pandas to match variants with annotation.
There are many databases. The biggest popular database is gnomad and the biggest disease database is clinvar. Others include omim, snpdb and 1k genomes. I'm not sure about blood sorry but I'm sure you can find one.


text:
correlation between imputed genotype and true genotype
There is a graph in the link https://odelaneau.github.io/GLIMPSE/rsquare_eur.html. This graph shows the imputation accuracy of certain MAF bins. Now I understand that $r^2$ is the correlation between imputed genotype and true genotype. What I don't understand is how they calculate this $r^2$ correlation. Can anyone suggest to me how I can define the r^2 correlation between true genotype and imputed genotype? I suppose it depends on the individual who defines his own $r^2$. So, it may differ for the different algorithms. However, different aspects or different $r^2$ definitions will help me understand and represent this graph slightly differently.

This is done by downsampling. Take the 1000 genomes, set some genotypes as missing ./., impute them using GLIMPSE, then measure correlation between the genotype dosages in the imputed dataset and the same sample at high coverage - in this case, the 1000 genomes sequenced to 30x coverage by the NY Genome Centre.

Is the correlation a Pearson correlation coefficient?

@winni2k judging by this yes I think so https://github.com/odelaneau/GLIMPSE/blob/master/concordance/src/utils/basic_stats.h#L110


text:
Tutorials on phasing and imputing low-coverage sequencing data
I am new to low-pass whole genome sequencing and have the basic idea of phasing and imputation. I have .vcf file after calling haplotypecaller tools from GATK. After searching the phasing and imputation tools I've encountered some tools like Beagle, Minimac, Shapeit, Glimpse, Eagle. The problem is that I didn't find enough hands-on material on how to process a .vcf from unphased to phased and then to an imputed one.
Can anyone suggests any hands-on tutorial regarding phasing and imputation?

What is your sample size, organism, and average coverage per sample?

total sequences read 4093666, Human cell, average coverage (0.1x-6x).

ok, here is the tricky bit: are you trying to impute germline or somatic mutations? I have tried beagle with somatic mutations, and it simply doesn't (and shouldn't) work. If you are using germline mutations, what is your reference panel size?

I am trying to impute germline. I'm trying to impute it with 1000kg phase3 reference panel.

1000 genome phase3 reference panel.

This is all relevant for data 0.5x-1x coverage.
Assuming you have genotype likelihoods data, if you want to phase low-coverage data, the most suitable options are GLIMPSE, Beagle4 and QUILT. Of the three, GLIMSPE and QUILT are the most recent and similar in performance. shapeit, eagle and minimac aren't suitable for low-coverage data.
I have found GLIMPSE to be fast and memory efficient, but it requires several different steps; breaking the genome into chunks, imputing, reforming the chunks and then sampling from the haplotypes to phase the data. It has a fairly good tutorial here.
QUILT has fewer stages and faster, but is unlikely to make much of a difference unless you are working on very large datasets. It has a tutorial here.
Which one you use probably doesn't matter a whole lot, but I would probably suggest using QUILT as it has fewer stages and is slightly more recent.

I am interested in using Beagle. But it doesn't have good documentation and tutorial. Can you suggest any?

Only beagle4 is suitable for low coverage data and it’s pretty old and doesn’t scale well, so I’d advise against using it unless you have a good reason to?

Also, I guess Shapit is only for phasing. I found some tutorial on Shapit+impute2 https://mathgen.stats.ox.ac.uk/impute/impute_v2.html#examples. What do you suggest about it?

Glimpse or Quilt which one is much older?

Also, there is a new Beagle5.2 version. Will it be suitable for low pass sequencing?

GLIMSE is slightly older, Beagle5.2 is not suitable for low coverage and neither is the shapeit2 and impute2 tutorial you linked.

Thanks for your comment. But I found a paper on the performance of phasing and imputation tools https://www.karger.com/Article/Fulltext/489758. The data here has pretty low depth(4x). It suggests that beagle 4 or shapeit+impute2 has a good performance. Can you tell me more elaborately why are you telling that shapeit+impute2 is not a good choice for low pass sequencing? Maybe add your suggestion in the answer?

I suppose it depends on what you consider 'low-pass'. I am referring to data around 0.5-1x coverage, so the results from that paper aren't really relevant to that case, and it predates the release of GLIMPSE and QUILT. Please read the papers I linked for more details. I think that’s everything I’ve got to say about this, as comments are not for extending questions and discussions.

That said, if you want to add in more specific detail about the coverages you are interested in, I will amend the question.

The author of shapeit (Olivier) once told me that shapeit2 imputation simply doesn't work below 4x. I have no reason to doubt that assessment.

I investigated a few different ways of running beagle v3 (same algorithm as v4) for this paper https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4522619/ i found that beagle scaled really poorly in run time with ref panel size. This is why we went just with asian 1kg haplotypes. I think beagle could work for you. There's also quite a bit about that testing in my thesis: https://ora.ox.ac.uk/objects/uuid:19ce0c8e-1d65-44ff-b56c-77ee849b2167

That said, the imputation using beagle on 11k whole genome samples (1.7x) took around 300k cpu hours. I would expect quilt to require less compute time and probably get you similar results.


text:
Interpreting imputation result from GLIMPSE
I'm following this tutorial of GLIMPSE for learning. I was expecting some extra SNPS coming from the 1000 genome reference in the resulting .vcf file. Though I understand the phasing in the output resulting .vcf, what I don't understand is the imputed result. Apparently the .vcf before and after has the same number of SNPs. I ran the scripts from here. I'm including the vcf files before and after imputation.
Can anyone explain me where the imputed SNPs are in the output v.cf file. If my concept is wrong can anyone tell me what should we expect in result.

Please share the exact command you used.

Just run the scripts from https://github.com/odelaneau/GLIMPSE/tree/master/tutorial

Please paste them into the question to make it easier for myself and others to know exactly what you did. Questions should be as self-contained as possible, a external links break.Thank you.

It seems that you didn't use the following flag:
 --impute-reference-only-variants   Allows imputation at variants only present
                                     in the reference panel. The use of this
                                     option is intended only to allow
                                     imputation at sporadic missing variants.
                                     If the number of missing variants is
                                     non-sporadic, please re-run the genotype
                                     likelihood computation at all reference
                                     variants and avoid using this option,
                                     since data from the reads should be used.
                                     A warning is thrown if reference-only
                                     variants are found.

You were expecting additional SNPs which were present in the reference panel to be imputed into your targets, but you didn't include the flag --impute-reference-only-variants, so only individual missing genotypes were imputed. Individual missing genotypes ./. are always imputed if they are present in both the target and reference, but if the SNP is present in the reference and not the target, it won't be unless you use that flag.

Doesn't answer my question. Can u explain it clearly?

You mean without --impute-reference-only-variants  it will give me only phased data of the sample? Not the imputed snps from reference?

Can you explain your answer more elaborately? I didn't fully understand the answer .

@ShafayetRahat With all due respect, the function of SE isn't to provide a tutorial service for everything you might not totally understand, especially when you say 'I don't understand' rather than being specific about what you don't understand. I feel like my answer provides a clear answer to the question you gave. Having to respond to multiple additional questions in the comments gets a bit tiresome after a while and doesn't;t motivate me to answer other questions you might ask.

sorry about that. Actually, I'm at learning period. So, it may seem that I was asking for a tutorial service. But, in reality, I'm just stuck in the process. I guess the purpose of the whole stack exchange community is for helping with a problem to solve it. It's not like you are bound to answer the question or something like that. You can choose not to answer the question. Anyway, my point is I'm not treating anyone as my personal instructor. I'm very sorry if any of my conversations hurt anyone. And once again sorry for your inconvenience.


text:
Where can I download 30x 1000 genomes cram files?
From the preprint published by 1000 genome project (https://www.biorxiv.org/content/10.1101/2021.02.06.430068v1.full) I think the 30x data is for WGS. Can anyone confirm for me if the following file is really WGS data?
ftp://ftp.sra.ebi.ac.uk/vol1/run/ERR324/ERR3242643/HG03007.final.cram.
Also, Is there any more freely available sequencing data where I can find WGS data for 30x and low coverage(<10x) for the same sample for the South Asian population?

Yes, those are the .cram files for the 30x whole genome sequencing data for 1000 genomes.
If you want 10x and 30x sequence data for the same individual, you can just downsample the 30x .cram file using samtools view -s.

30x WGS or WES? Can u specify it in your answer?

WGS - I don't know why you'd think it was exome sequencing - there's no mention of that in the paper or the website.

There is actually. But that was at the depth of 65.7x coverage. However, I asked just to be sure, anyway once again thanks for your answer.


text:
How to use combat in order to remove batch effects?
I have RNA seq data and I need to use combat to remove the batch effects. Somehow when I run it, it isnt actually doing anything.
The code:
dds <- DESeq(dds, betaPrior=TRUE)
normalized_counts <- counts(dds, normalized=TRUE)
log2 = log2(normalized_counts+1)

modcombat = model.matrix(~dex, metadata)# - with metadata being a variable containing treatment or control (under dex column) and batch and name of each patient.

com<-ComBat(log2, metadata$Batch, mod = modcombat)

Its supposed to be 4 different batches, but in the com variable I can see that the values have stayed the same as in log2.
What could be wrong? Would appreciate any help!

side note: combat is a linear batch correction method and does not perform optimally on all datasets

Are you completely sure that you can make the model matrix with a mix of column names as strings and object names like that?

No I am not sure about that, I havent found a good explanation for how to make the model matrix. I suspect my problem is there maybe.


text:
How do I install hammock?
I'm working with phage peptides libraries. Currently I'm analysing a massive data set from phage selection, NGS readout, preprocessed already (fasta file). I want to analyse the sequences to cluster peptides according to their consensus motifs. I found that Hammock may be helpful to do this.
Here is the GitHub repo for Hammock: https://github.com/krejciadam/hammock. I downloaded the master, unpacked it, but I don't know how to use it.
I guess I need to compile, but I have no idea what files need to be compiled, how to do it, what software I need use to compile etc.
I work with Mac/terminal.

Download a release (Hammock_v_1.2.0.7z) from https://github.com/krejciadam/hammock/releases
Unzip
To run the first example:

cd /your_path/Hammock_v_1.2.0/examples/MUSI/
java -jar ../../dist/Hammock.jar full -i musi.fa

As mentioned in the comments, it's in the manual and any further issues need more information (errors etc).

Answer from @ram-rs, converted from comment:
There's a Hammock.jar file - there' no need to compile Hammock, although they recommend compiling dependencies. Have you tried running java -jar dist/Hammock.jar --help? Follow their steps from Example 1 - it should work.


text:
Find corresponding symbol for gene used in Seurat
I have a Seurat object
> pbmc
An object of class Seurat
48608 features across 11709 samples within 2 assays
Active assay: SCT (22220 features, 3000 variable features)
 1 other assay present: RNA
 2 dimensional reductions calculated: pca, umap
>

I want to plot TCRVδ2
> FeaturePlot(pbmc,"TCRVδ2")
Error: None of the requested features were found: TCRVδ2 in slot data
In addition: Warning message:
In FetchData(object = object, vars = c(dims, "ident", features),  :
  The following requested variables were not found: TCRVδ2
>

This data base is about human TCR genes
https://www.genenames.org/data/genegroup/#!/group/370
But I even can not find this gene here
Do you know what is the corresponding symbol for this gene used in Seurat?

TCRVδ2 is not a gene, but it defines a specific subset of gamma-delta TCRs. You can look at TRDC (delta chain) and TRGC1 and TRGC2 (gamma chains). TCRVδ2 should have the combo TRDC-TRGC1, while TCRVδ1 should express TRDC-TRGC2. You may need to increase the resolution of your clustering, to avoid clustering the gamma-delta T cells with the NK cells.

Thank you so much. I have another question about gene symbols. CD45 is corresponding to PTPRC butler most of CD genes I can not find the corresponding official gene symbols. Please, do you know how I can find corresponding gene symbols for CD genes?


text:
How do I install shapeit4 on Ubuntu/Mint?
I'm having trouble install shapeit4 4.2.2 on my new Mint 21 installation (Ubuntu 22.04 LTS). According to the official docs, the requirements are:

HTSlib: A great C library for reading/writing high-throughput sequencing data.
BOOST: A free peer-reviewed portable C++ source libraries. SHAPEIT4 uses two specific BOOST libraries: iostreams and program_options.

I installed HTSlib with make and sudo make install, and libboost using its calls ./bootstrap and ./b2 but I still get errors when trying to run the make file for shapeit4:
/usr/bin/ld: cannot find /home/user/Tools/htslib-1.11/libhts.a: No such file or directory
/usr/bin/ld: cannot find /usr/lib/x86_64-linux-gnu/libboost_iostreams.a: No such file or directory
/usr/bin/ld: cannot find /usr/lib/x86_64-linux-gnu/libboost_program_options.a: No such file or directory

Manually editing makefiles is just... ugh. Save yourself the effort and use conda:
conda create -n shapeit
conda install -c bioconda shapeit4 -n shapeit

The documentation appears to be missing an edit one has to do. Open the makefile in the shapeit4 dir, and you will see something like (for version 4.2.2):
#COMPILER MODE C++11
CXX=g++ -std=c++11

#HTSLIB LIBRARY [SPECIFY YOUR OWN PATHS]
HTSLIB_INC=$(HOME)/Tools/htslib-1.11
HTSLIB_LIB=$(HOME)/Tools/htslib-1.11/libhts.a

#BOOST IOSTREAM & PROGRAM_OPTION LIBRARIES [SPECIFY YOUR OWN PATHS]
BOOST_INC=/usr/include
BOOST_LIB_IO=/usr/lib/x86_64-linux-gnu/libboost_iostreams.a
BOOST_LIB_PO=/usr/lib/x86_64-linux-gnu/libboost_program_options.a

So the paths appear to be hardcoded. Set these to the paths where you installed htslib (and compile it first).
For libboost, I found that its own installer doesn't place files in this expected location. So you can either change these paths to where you installed libboost, or you can use the apt manager version which will put the files in the right location by default. This is done with sudo apt-get install libboost-all-dev (found here).
After this, just run the make file again, and the output should looks like this:
user@computer:/dir/shapeit4-4.2.2$ make
make: Nothing to be done for 'all'.

Actually, precompiled SHAPEIT 4.2.2 is available in a Ubuntu repository.
sudo apt install shapeit4

Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.


text:
Time for running ADMIXTURE analysis
I am trying to run the ADMIXTURE software for the first time to analyse the structure of an in-house dataset of samples comparing to 1000 genomes project ancestral populations. I am trying to determine ancestry just from the X chromosome of those individuals.
To do that, I have downloaded the X chromosome file from 1000 genomes (ALL.chrX.phase3_shapeit2_mvncall_integrated_v1b.20130502.genotypes.vcf.gz). I have removed duplicated with bcftools and converted to plink binary with the following command in plink 1.9:
plink --bcf ALL.chrX.phase3_shapeit2_mvncall_integrated_v1b.20130502.genotypes.bcf --keep-allele-order --vcf-idspace-to _ --const-fid --allow-extra-chr 0 --split-x b37 no-fail --make-bed

Next, I modified the .bim files to change the polymorphisms ID so them where the same as their position, as a way to standarize their names. Then I filtered related individuals, also with plink and found the intersection of .bim files from 1000 genomes and my in-house dataset (a bed file with microarray genotyped variants only from X non PAR region), that was already filtered to contain only variants genotyped against the plus strand (without indels).
Then I merged both datasets with the following plink command:
plink --bfile filename1 --bmerge filename2.bed filename2.bim filename2.fam --make-bed --out filemerged

That command generated a file with all the SNPS that would have more than 3 alleles in the merged file so I removed all them from both datasets and merged them again. On that merged file I pruned SNPs by linkage disequilibrium:
plink --bfile cromXmerged2 --indep-pairwise 50 10 0.1

, assign sample sex on the .fam file and filtered individual with missingnes higher than 10%:
plink --file mergedfile --mind 0.1 --make-bed --out

I also set all heterozigous calls as missing with --set-hh-missing plink command.
With the final input I runned a supervised ADMIXTURE as following:
./admixture finalinput.bed 3 --haploid="male:23" --supervised -j4

Before running of course I created the finalinput.pop file in open office calc, by creating a list with the assigned population for each sample, with my in-house dataset left blanc. I just used african, european and east asian 1000 genomes samples because my own dataset has brazilian admixed samples. The final dataset has approximately 1500 SNPs and about 1600 individuals. My computer set up is CPU: i3-2370M with two cores, 4 threads@2.4 GHz, RAM 8 GB@1333 MHz and an SSD 120 GB Sata 3.0. OS Ubuntu 20.04 64 Bits.
The thing is that ADMIXTURE have been running for three days without showing any error or producing result yet. I would like to know if it is normal to this analysis to take so long or if it could be any problem with my dataset.
I would appreciate any suggestion or advice!!!
Thanks in advance,
Iriel

This should take a few minutes, not 3 days! Definitely something seems to have gone wrong here. Are you running it locally or on a remote server? ADMIXTURE produces an output as its running, so you should check that it's doing that

I am running it locally. And it did not produced any output at the moment. I tried a unsupervised analysis before with K = 3 to 10 and it took a few days to produce Q and P outputs for the first 2 K values. But I canceled that run because supervised analysis would be better for my goals.

Did it literally not produce any output at all? I have just had a quick check and run it, and it produces this almost immediately:
Parallel execution requested.  Will use 8 threads. Random seed: 43 Point estimation method: Block relaxation algorithm

Nothing, it keeps showing these lines:743116 (QN/Block)  Elapsed: 0.354 Loglikelihood: -1.1648e+06 (delta): -103.745
743117 (QN/Block)  Elapsed: 0.331 Loglikelihood: -1.16479e+06 (delta): 6.94675
743118 (QN/Block)  Elapsed: 0.327 Loglikelihood: -1.16472e+06 (delta): 67.0046

I mean, It did not produce any file yet.

Seems like it's having trouble converging. There is a post here (https://www.biostars.org/p/368113/) which suggests that using  haploid mode and multithreading might cause issues converging and that you should try running it without multithreading. Actually, your dataset is so small that I doubt you need to multithread anyway :)

Thanks for the reply! I will try without multithreading then.

no problem. I am interested to see if this solves the problem.

Yes! I took just 6 seconds to run!

Great. Should add the solution as a response to this and accept the answer so other people who had a similar problem can find the solution.

Apparently there is a bug in ADMIXTURE that prevents for converging when using at the same time the haploid mode and multithreading (flags --haploid and -j respectively). The problem is solved by avoiding using one of them.


text:
NO_COOR reads not in a single block at the end 0 -1
I am doing an RNA seq analysis with mouse samples. I had my reads in several lanes, so I concatenated those fastq reads. Then I've mapped my reads with Hisat2 against the new version of mouse  "GCA_000001635.9_GRCm39_genomic.fna".
Apparently everything works well up to here. Now I'm using samtools 1.9 to convert sam to bam, to sort those bam and create the index, just like this:
samtools view -bS ${file} > ${file}.bam

samtools sort -n ${file} -o ${file}_sorted.bam

samtools index  ${file}_sorted.bam

But when indexing I get this error:

[bam_sort_core] merging from 11 files and 1 in-memory blocks...
[E::hts_idx_push] NO_COOR reads not in a single block at the end 0 -1
samtools index: failed to create index for
"CD4_1A_concat.bam_sorted.bam": No such file or directory
[bam_sort_core] merging from 11 files and 1 in-memory blocks...
[E::hts_idx_push] NO_COOR reads not in a single block at the end 23 -1
samtools index: failed to create index for
"CD4_1B_concat.bam_sorted.bam": No such file or directory
[bam_sort_core] merging from 11 files and 1 in-memory blocks...

And more...
Any help would be appreciated!

This is essentially a duplicate of Samtools Index: Chromosome Blocks not Continuous but with different error messages, and a duplicate of Error given while trying to index a BAM file with Samtools Index - NO COOR?.
samtools sort -n ${file} -o ${file}_sorted.bam

These error messages also indicate that the reads are not sorted by coordinate — in particular, that the unplaced unmapped reads with RNAME * are not all together at the end of the file.
This is because samtools sort -n has been used to sort the reads by name instead. Remove -n to sort by position, which is what is needed to prepare a BAM file for indexing with samtools index.

Thank you! problem solved :)


text:
About the log2 fold change
It seems that we have two calculations of log fold change:
Actual log2(FC) = log2(mean(Group1)/mean(Group2))
Limma's "Log(FC)" = mean(log2(Group1)) - mean(log2(Group2))
For volcano plots, can we use the second one?
https://www.biostars.org/p/100460/

Aren't the outputs of both calculations the same?

log(x/y) = log(x) - log(y) -> this is log math. Like @RezaRezaei says, the two calculations are the same. I guess there could be differences owing to how computers calculate the values.

@RamRS Could there be a small difference due to finite precision of log transformation?

@KamilSJaron I guess it could be a factor. I don't know how computers do complex math, but I'm guessing that's where the difference lies. I think for x=20, y = 5, the two log2 values should be the same. I do think that even for more complex math, precision should be maintained to the third decimal place at least, so how significant would the difference be?

Yes, you can use the second one for volcano plots, but it might help to understand what it's implying.
The difference between these formulas is in the mean calculation. The following equations are identical:
log(a/b)

and
log(a) - log(b)

However, in your example, the first equation calculates the mean prior to doing the log, whereas the second example does the log first (i.e. log(mean(Group1)) vs mean(log(Group1))).
As a consequence, the equations are similar, but not identical, and that difference depends on the underlying assumptions of the expression model.
In the first equation (what you call 'Actual'), the genes are assumed to follow a normal distribution in linear space, whereas in the second equation they're assumed to follow a normal distribution in log space.
Which one of these is correct is a subject of debate. It's pretty much established that gene expression for all genes together is log-normal (or, at least, much closer to log-normal than to linear), but whether this applies at an individual gene level is less clear.
In most cases, the results will be similar when looking at individual genes, and the equations can be used interchangeably. Where they differ substantially, there will be large expression differences between genes, so it's unlikely to alter the interpretation of results.

You can use what limma reports for everything, including volcano plots. In the post you link to I describe in a comment why this is done in every package that I know.


text:
R not being able to find the file despite being in correct directory
I have question regarding a project I have to do about a health student study. I cannot import the csv file via read.csv or read.csv2 despite having transformed all files to csv and being in the correct folder.
Here are the command lines:
setwd ("C:/Users/fek/Desktop/Analyse i-Share")
getwd()
TA<-read.csv2(file=Biobanque.csv, header=TRUE, sep=",", dec=",", fill=TRUE)

And here are the results in the console window:
setwd ("C:/Users/fek/Desktop/Analyse i-Share")

getwd()

[1] "C:/Users/fek/Desktop/Analyse i-Share"
TA<-read.csv2(file=Biobanque.csv, header=TRUE, sep=",", dec=",", fill=TRUE)

[1] Error in read.table(file = file, header = header, sep = sep, quote = quote,  :
objet 'Biobanque.csv' introuvable
To translate from french, it says it can't find the file!

The problem is just what it says.  By omitting the quotes around "Biobanque.csv", R thinks you are talking about an object whose name is Biobanque.csv. And R is confused because of course you have not created an object by that name.

This is correct - the error message states that R cannot find the object 'Biobanque.csv', not that it cannot find the file 'Biobanque.csv' (which would give the error "No such file or directory").


text:
How can I subset WGS data to the level of WES variants?
I would like to compare mutational signatures1 in patients from different studies, however some studies are based on exome seq (i.e. ~20,000 coding variants) and some are from whole genome seq (i.e. ~22,000 coding variants) - is there a way to 'downsample' WGS data to better reflect WES data and effectively 'ignore' the coordinates of those ~2000 coding variants in the VCF files?

1Alexandrov, L.B., Kim, J., Haradhvala, N.J., Huang, M.N.,
Ng, A.W.T., Wu, Y., Boot, A., Covington, K.R., Gordenin, D.A.,
Bergstrom, E.N. and Islam, S.A., 2020. The repertoire of mutational
signatures in human cancer. Nature, 578(7793), pp.94-101.

Apparently this is trivial - according to a colleague you can take the WES bed file coordinates and remove any SNVs in the WGS file that fall outside those regions - but he didn't know of any software that does this in a more sophisticated manner.

bcftools would be my choice, I'm sure bedtools could do the trick, too.
Something along this
bcftools view --regions-file

or  --targets, --regions-file might require a tabix index.

I'll expand slightly on the previous answer. First print off the positions for the exome file and then use bcftools view to filter the variants from the whole genome file. You could also index the whole_genome.vcf file to make the filtering faster.
bcftools query \
    -f'%CHROM\t%POS\n' \
    exome.vcf > exome_variants.txt

bcftools view \
    -T exome_variants.txt \
    whole_genome.vcf > whole_genome.exome_positions.vcf


text:
Denovo, Stacks: Getting an “ambiguous redirect” error
I tried running the following command
/public/home/WHC_BBT/Softwares/bin/bin/denovo_map.pl --samples $RawData -O $Popmap -o $fullM2n3 -M 2 -n 3 -m 3 -T 15 &> $logfile

This is the error I got

/opt/gridview//pbs/dispatcher/mom_priv/jobs/14923.admin.SC: line 6:
$logfile: ambiguous redirect

I need help on how to resolve this.
Thanks.

This most likely means your variable $logfile does not exist.

At a guess, the variable $logfile is not properly defined (i.e. it hasn't been defined at all). What is the output of the following command:
echo "logfile: ${logfile}"


text:
How to select specific lines and specific columns
I have  a data frame that contain lot of columns and lines, like this for example
Column1   Column2
5            T1
1            T2
1            T3
1            T1
2            T1
1            T1
2            T1
1            T1

I want to extract just lines that have the value 1 in column 1 and Value T1 in Column 2
So the result I want is
Column1   Column2
1            T1
1            T1
1            T1

Can you tell me please how I can do that

Hi @James what code have you written so far? You must have imported into R and done some basic dataframe work.

Hello , thank you very much for you answer , yes i tried this one : Pokemon_2 <- dropRows(ds, Pokemon$Generation == 1 , Pokemon$Type_1=="Water" )   but it didn't work  , Error in dropRows(ds, Pokemon$Generation == 1, Pokemon$Type_1 == "Water") :
could not find function "dropRows"

You should look at the subset function. Please google "filter data frame by column value" - much like writing for loops in basic programming languages, it is a basic operation that does not warrant creating new posts on forums.

I’m voting to close this question because it is a basic R operation has nothing to do with bioinformatics specifically

You can simply use subset function in R and use the code given below,
    new_dataframe <-
    subset(name_of_your_dataframe,Column1 == 1 &&
    Column2 ==  "T1")

Please keep in mind that R is case sensitive and you should be careful while using the variable names.


text:
Loss of predictive power of polygenic risk score when dataset contains missing variants
I am trying to calculate polygenic risk scores (PRS) scores for a new dataset. This dataset does not have all the variants that the PRS score needs. The PRS score I am interested in has 40 variants, but in my dataset only 20 of those variants are present. What is the correct way of assessing the decreased predictive power of the PRS score due to the missing variants?
I assume it is not as easy as summing the effective weights of those variants I have and comparing to the full sum (all variants needed for the PRS).
Thanks a lot and sorry for newbie question!

In a sample with all 40, determine the accuracy of the PRS. Then in the same sample, remove the 20 you are missing in the other set, and then determine the accuracy with the missing variants. This will give you an idea of the level of attenuation. However, you may not have access to such a dataset. In that case, I would recommend imputing your data before calculating a polygenic score. (It's likely that the GWAS the PRS is from used imputation itself.)

Which software package around you using, PRSice? Generally missing values comprising half the total observations are extremely difficult to calculate using regression analysis (or many statistical analyses for that matter.
One way forward is to include (an) example(s) of all 40 variants into your data set and perform the analysis with and without the added data, then compare this against a data set where all 40 variants are observed and your samples are excluded. You might find the missing values have low probability of PRS and this would then explain your results.
Missing analysis in regression is a specialist area of statistics and is not for the faint hearted, so a work-around is the best way forward.


text:
Inheritance annotation for nsSNPs
I am relatively new to SNP analysis. Is there a database to find annotation about the Mendelian inheritance of SNPs?
I have a small list of nsSNPs (non-synonymous SNPs) and I need to find how they are they passed to from parents to offspring.
Thank you in advance.

If you only have a small number of genes, OMIM is a very easy way to do this.
For example:

If you have too many to do manually, I suspect there is an easy way to query lists, although I'm not sure exactly the way myself. Perhaps someone can add a new answer if so.


text:
Single-cell RNA sequencing proper abbreviation
I wonder which form of a single-cell RNA sequencing abbreviation is correct:

scRNA-seq

scRNA-Seq

In the same single-cell sequencing Wikipedia article, both forms appear.
There is also inconsistency among papers I've seen - some authors prefer the first form, others the second one.
If there's any reason behind the different form, I'd like to hear about it.


text:
16S merging sequences, removes ALL reads
I am currently analysing some 16S (V4 region) rRNA data, using the qiime2 command-line software. Everything is working fine for the most part, until we are to merge the forward and reverse sequences - this is where the problem arises.
The command used to do this, is the following:
qiime vsearch merge-pairs \
  --i-demultiplexed-seqs reads_qza/reads_trimmed.qza \
  --output-dir reads_qza/reads_joined

Where reads_trimmed.qza, is a qiime2 'artifact' i.e. object, which contains all of our samples which have been run through cutadapt, in order to trim primers/adaptors.
As the title says, when merging the pair of sequences, we get almost no sequences merged, in the end, please refer to this image for a visual representation of what this looks like:
We tried to modify the above script to be more lenient with regards to the merging of sequences, by adding the argument --p-maxdiffs 50, which is defined as:
Maximum number of mismatches in the area of overlap during merging.
But to no avail, we still have basically no matches, with a median read number of less than 1.
What could be causing this issue? Has anyone faced this before? Could this have to do with the sequencing technology? We used Illumina NextSeq, and because of it opted to use qiime and deblur rather than dada2, due to binned quality scores associated with this technology, there's an interesting discussion on this here:
https://github.com/benjjneb/dada2/issues/1307
Thank you,

You're attempting pair-ends assembly without complete coverage, e.g. short-read Illumina.
Switch to fastp for QC and follow this up with Megahit or metaSpades. You will lose sequences maybe by as much as 10-fold loss, but you'll recover quality data.
If this fails there are alternative approaches.

Thanks for your response. My understanding was that these methods (qiime2) were designed for short-reads though, were they specifically designed for long reads ? I will look into these other methods, though !

qlime2 claims are clearly not correct.

Noted ! Thanks for your input.


text:
Co-occurrence networks in Metagenomics studies
I have recently acquired some 16S metagenomics data, and was wondering if anyone can speak of the potential limitations, challenges as well as advantages to conducting a network-based study on metagenomics, such as what is done w.r.t co-occurrence networks.
I thought it would be interesting to deviate from the "classic" study of taxonomic abundancies present in different timepoints / conditions. However, it seems as though for such studies they are usually dealing with a relatively large number of samples - please find paper below:
https://animalmicrobiome.biomedcentral.com/articles/10.1186/s42523-020-00054-4
Could such analyses be done with much smaller sample sizes such as those generated via mouse models? What are the potential limitations of such an approach, and how could they be overcome?
Thanks.

Do you mean covariance, ie this is the measurement of co-occurance?

No. There are many possible ways to generate the networks. Based on my understanding, they are generated primarily using correlation metrics, in order to capture 'pairwise affinities', between different taxonomic classifications - i.e. bacteria i, and bacteria j.

I do agree, its just 'co-occurance' isn't a statistical method. Pearson's or Spearman's correlations for example? Pairwise is a statistical method but not a metric and might even suggest distance method like P-distances :-/. Without knowing the metric you run the risk that the assumptions for one, may not be applicable. At present its not clear whether the analysis is applied directly to sequence data - which would likely be a mistake - or their taxonomic designation.

In one of the studies I am currently reading, they use a 'treshold-Kendall correlation'. My understanding is that this is to be applied to processed sequencing data, after taxonomic classification, rather than simply the sequences.

Okay, that makes sense. Would be good to update your question. Its Spearman type stuff which is part of SciPy stats. I dunno enough about R, but there are some very capable R-coders here. Really this is population dyanamics of a microbial community, lactobaccilus (if its gut) and all that. Its definitely not my thing and I wouldn't approach the question in this way.

Would you care expanding on what you don't like about this particular approach? Doesn't have to be a long answer, I just want to get an idea of the limitations of this approach ..

To answer your question more formally. I strongly believe you are treating the mice with something that will shift its microbial flora. Clearly the experiment is not described.
Thus what I simply feel is you essential have a discrete time series comprising before, during and after treatment. 'Treatment' could be anything drug, vaccine, antibiotic, diet who knows.
I think that is a very different time-series to a one where classic correlation statistics will be informative, I'm just guessing in the absence of a hypothesis, thus I could be wrong.
I also kinda suspect lactobacillus are strong players.
Anyway the stat you are looking at is a special case of the general correlation coefficient in much the same way Spearman's classic statistic is. What I can say is that the statistic is not a specialist approach for large sample sizes, which I believe is your question.
'New fangled cool', might be an approach of data mining and that might be okay if you've loads of data you don't know how to analyse, like it goes this way, that way. Ultimately, statistics are used to assess an exact hypothesis. If e.g. lacks power (very common) or the assumptions are violated then look deeper. I do agree that in population biology generating a hypothesis is often helped by data mining, but what I think is you already have an exact hypothesis and I'm far from sure this would nail it.


text:
Frequency of specific viral sequence in .BAM or .fastq
I was wondering if anyone had any experience 'counting' the frequency of a particular viral sequence in an individual's fastq sequence. So basically, counting the number of occurrences of a particular sequence given someone's fastq file.
The thing is, we must speak with the sequencing company regarding what format we'd like the files and must choose between .bam and .fastq - both would be extra $. Usually I'd choose .bam out of convenience but I was wondering, given that it has already been aligned the human genome I might not be able to proceed with this 'counting' process. So, my question is: if I were to use samtools - or any other program - to convert the .bam files to .fastq, is there any possibility that I could lose such information - e.g. viral seqs - that would have otherwise been available should I have just gotten the .fastq files directly?

Rule of thumb is: never go with FASTQ if you can avoid it, it’s a net waste of space. Post-alignment BAMs are as small, or smaller, than gzipped FASTQ, despite containing more information.

Most read aligners will report unaligned reads as well, which presumably will include your viral sequences. I would ask them to formally confirm that the BAM files will contain unaligned reads before choosing that option.

Thank you !! Will check with the company to be sure

If you are interested in a particular viral sequence, it seems wasteful to request bams that are aligned to the human genome, but you could request bams aligned to the human genome with unaligned reads retained. As an alternative, you can generate unaligned bams, where all of the same reads are present that would be present in a fastq but in a bam format. You can revert this to a fastq or align the fastq directly to a viral reference yourself. Aligning to small references is fast.


text:
Find common SNPs in multiple VCF files
I have 3 VCF files.
A1.vcf
A2.vcf
A3.vcf

I want to get the common SNPs that are present in all these three files. And output must be in vcf format.
output:
common_A.vcf

Do you mean SNPs in terms of shared variants or shared positions?

google intersect or overlap between VCF files, -> bcftools isec

I have tried isec, but the produced files do not contains the VCF file contains common file

A potential solution is to compress and index the files using bgzip and tabix, then use bcftools isec to find SNPs common to all three files, e.g.
for files in A1.vcf A2.vcf A3.vcf
  do
  bgzip "$files"
  tabix "$files"
done

bcftools isec -n=3 -c all -o all_common_SNPs.vcf A1.vcf.gz A2.vcf.gz A3.vcf.gz
```

Hi @Jared_mamrot - can you please explain -c all I've seen -c none -c all/both I'm unable to understand/follow this flag.

As I understand it @Death-Metal, if you have two vcfs, one with the variant T->CCA at chr1 pos10001 and the other with the variant T->A at chr1 pos10001 , the -c all flag will treat these as a 'common to the two files' (i.e. there is a variant at chr1:10001 and it will be 'kept' when you isec). If you specify -c none, the REF/ALTs need to match exactly (i.e. both vcfs need to be e.g. T->A to be kept). With -c both, the variant must be the same type e.g. both SNPs or both INDELs, regardless of the ALTs. More details here: https://samtools.github.io/bcftools/bcftools.html#common_options


text:
Trim Seurat object to contain expression info only for selected genes
I'm working with some large Seurat objects (MOCA, MCA, Tabula Muris) studying gene coexpression, and I'm running into memory issues.
Is it possible to remove all genes in a Seurat object that are not in a specified gene array? Ideally I'd like to be working with a small R matrix with say, 250k rows (cells) and 10-50 columns (gene expression data) which will enable fast processing and comprehensive analysis using other packages.

There might be drawbakcs with this approach, for example if you end up with genes that are highly expressed in one cell type and lowly expressed in the other, these two sets of cells will look like expressing these genes at similar levels after normalization if normalization is performed with this limited number of specific genes. You might want to consider sampling your cells, you can decrease the number of cells instead of genes.

Is your memory usage low enough that you can initially read everything in and then filter later? If so you could read the Seurat object in like normal and then use subset against a vector of your gene names.
subsample <- subset(seurat_obj, features = my_genes)

A different approach if you are using Seurat3, is DietSeurat(). It allows you to diet the object by removing the components that you don't need. For example you can keep the normalised/scaled matrix and remove the raw counts. This approach could reduce space and memory usage, while keeping all your genes in place.


text:
Seurat VlnPlot presenting expression of multiple genes in a single cluster
Seurat VlnPlots are most commonly used to visualize differences in any given gene expression across multiple clusters or cell types. For example:
VlnPlot(object = MouseCellAtlas, idents = c("T cell","Neutrophil","Erythroblast","Monocyte","Macrophage"), features = 'Fzd9')

But suppose I want to view the expression of several genes in just idents = "T cell". I can do this to generate three separate plots:
VlnPlot(object = mca, idents = 'T cell', features = c('Fzd9','Ctnnb1','Apc'), combine = TRUE, ncol = 3)

But ideally I'd like this second plot to look just like the first plot, except with only a single tissue on the legend ("T cell") and multiple genes in the various colored violins. Unfortunately, I am not aware of any simple way to combine these violins into a single plot area.
This doesn't do what I need:
VlnPlot(object = mca, idents = 'T cell', features = c('Fzd9','Ctnnb1','Apc'), combine = TRUE, ncol = 3)

Is there built-in functionality in Seurat for generating a VlnPlot for multiple genes in a single tissue in a single plot? Or can anyone suggest a workaround?

To do so one workaround it to have your data in "long format" and then use the column that holds the "gene names" as the x variable while plotting.
You can use FetchData() to extract data from a Seurat object. VlnPlot's default is the data slot (of the active assay if using Seurat v3 I suppose). And you can specify which cells and genes to retrieve.
selected_cells <- names(panc8$celltype[panc8$celltype == "gamma"])
data <- FetchData(panc8,
                    vars = c('FZD9','CTNNB1','APC'),
                    cells = selected_cells ,
                    slot = "data")
> head(data)
           FZD9   CTNNB1      APC
D101_5        0 0.000000 0.000000
D101_43       0 2.007853 1.001958
D101_93       0 0.000000 0.000000

Melt() would transform your data into "long format". By not specifying any arguments, we push for all of the info in the three variables to be gathered in two columns:
long_data <- melt(data)
No id variables; using all as measure variables
> head(long_data)
  variable value
1     FZD9     0
2     FZD9     0
3     FZD9     0

> tail(long_data)
     variable value
1870      APC     0
1871      APC     0
1872      APC     0

ggplot2 is used to add "violin" and "jitter" layers. You can customize the output to look (exactly) like the VlnPlot() output.
ggplot(long_data,
       aes(x = variable, y = value)) +
  geom_violin() +
  geom_jitter(size = 0.1)

And here is the same graph generated by VlnPlot():

There are no "violins" as the counts are almost entirely zeros. And see how different the ranges are in the y axes of the VlnPlot. So the "tweak" I have presented here would only work for genes that are expressed at similar levels / similar ranges.

thanks, this works well! The general approach will be useful moving forward with a number of other analyses as well. Thanks!


text:
What should be the numbers in the fastSTRUCTURE file?
The structure file for fastSTRUCTURE seem to be a table with rows and columns but the cells are all numbers, while my genotypes are simply bases. How could I know what I should insert there?

What format are your genotypes in at the moment? There is likely some tool someone has written to convert it to the correct format.

Just table of 'A/C' bases

'Table' isn't a format. What are the rows? What are the columns? Is it genotypes or alleles - i.e. is there one or two rows per marker? the more info you give, the easier it will be for someone to help you. Could you add a sample of your genotype format to the main answer so others can see it.

Numbers, usually in the form of 0/1/2, is a common way many formats represent bases.
For example, if a SNPs has bases which are A/T, then these are encoded as 0/1, with the reference allele/base taking the 0 and the alternate allele/base taking 1.
However, if you read the documentation, it says:

Specifically, rows in the data file correspond to samples, with two
rows per sample (note that only diploids are handled by this
software), and columns correspond to SNPs. The first 6 columns of the
file will be ignored; these typically would include IDs, metadata,
etc. This software only handles bi-allelic loci. The two alleles at
each locus can be encoded as desired; however, missing data should be
encoded as `-9'.

So you should be able to use your base letters here and it will work OK as long as the columns and rows are formatted correctly.


text:
Subsetting out prefix from genes and merging
I'm working with PTX samples (human and mouse), they're aligned to both transcriptomes and so gene names contain hg38- and mm10- at the beginning of each gene. I have eliminated teh prefixes because I want to integrate this dataset with a dataset that does not contain prefixes. I have annotated each cell type accordingly and have subsetted out the mouse cells.
human <- grep(pattern = "^hg38-", x = rownames(x = scData), value = TRUE)
mouse <- grep(pattern = "^mm10-", x = rownames(x = scData), value = TRUE)

mouse_1 <- substring(mouse,6)

human_1 <- substring(human,6)

human_only <- subset(PTX_no_prefix, features = human_1)

mouse_only <- subset(PTX_no_prefix, features = mouse_1)

Idents(mouse_only) <- "cell_ann"
TME_PTX <- subset(mouse_only, idents = c("Fibroblast_mouse", "Myeloid_mouse", "Endothelial_mouse", "Pericytes_mouse", "Neutrophil_mouse"))

I'm now trying to merge two datasets together but receive an error.
combined_het <- merge(TME_het_PTX, y = skin_only_ann)

scData <- combined_het

# remove cells with <200 RNA molecules, or >6000 molecules, or >30% mito
scData <- subset(scData, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt_mouse < 30)

#scData <- NormalizeData(scData, normalization.method = "LogNormalize", scale.factor = 10000)
#scData <- FindVariableFeatures(scData, selection.method = "vst", nfeatures = 2000)

all.genes <- rownames(scData)
GLOBALS_MAXSIZE <- 4000
options(future.globals.maxSize = GLOBALS_MAXSIZE * 1024^2)
scData <- SCTransform(scData, vars.to.regress = "nCount_RNA", verbose = FALSE)
#scData <- ScaleData(scData, features = VariableFeatures(object = scData), vars.to.regress = c("nCount_RNA"))
scData <- RunPCA(scData, features = VariableFeatures(object = scData))

combined_het <- merge(TME_PTX, y = skin_only_ann)
Error in intI(i, n = x@Dim[1], dn[[1]], give.dn = FALSE) :
  invalid character indexing

I checked the gene names from each dataset to make sure they're in the same format and they are.
head(rownames(skin_only_ann))
[1] "Xkr4"    "Gm37180" "Gm37363" "Gm37686" "Gm1992"  "Gm37329"
head(rownames(TME_PTX))
[1] "Xkr4"    "Gm37180" "Gm37363" "Gm37686" "Gm38148" "Gm19938"

possibly related? https://github.com/satijalab/seurat/issues/4452


text:
Why is there gene without annotated location?
I came across this gene Pepc (peptidase C)
NCBI gene database has a page for it here https://www.ncbi.nlm.nih.gov/gene/109616
But it does not contain any information on its coding sequence or genomic location. It just says it sits on chromosome 1.
Can someone familiar with how the reference genome is annotated explain why this is?

I agree that this is an odd gene: no transcript, but it has mouse strains and supporting literature.

I like the question a lot. You can't have this sitting on arguably (sure EBI would disagree) the worlds most comprehensive molecular database, not as comprehensive as it thinks it is. This is one for escalation with NCBI and the authors. If this was EBML Emily would step in. Definitely upvoted, I would recommend site members upvote to help the OP recover their bounty costs. I am aware this sort of thing goes on ... I've encountered in a completely different area.

Let me know if you intend to escalate it, because I will otherwise.

@M__ Thanks for the comment. I'm not sure what you mean by escalating. Do you mean reaching out directly to NCBI?

Yes, you can approach the authors but it never results in success (I've tried - not on this example).

This gene/protein appears to have been described exclusively as an isozyme / electrophoretic mobility assay of a protein.
See for example this PepC mutant page:

The Pepc locus controls electrophoretic mobility of the enzyme peptidase C of red blood cells. This LT/J congenic strain carries the Pepca allele from C57BL/6 that determines a slow anodally migrating electrophoretic band.

The relevant publications describing mutants are from the 70s and 80s, e.g. here, when there were no genomic resources for mice.
However, linkage mapping of electrophoretic variants could be used to establish linkage to loci on chromosome 1, apparently.
Very likely PepC is encoded by some annotated gene on chr1, but no one knows which one exactly because they don't really know much about the protein other than its electrophoretic variability. But that's just a guess by me.
When the databases were set up, they had to accommodate 1) mouse mutant strains, and 2) published articles referencing the protein, so they did their best.
However, no one ever followed up to link everything together, so we just have these orphan entries.


text:
What is the origin of HIV1?
What is the origin of HIV1?
This is the alignment of a small part of the POL protein of the HIV1 virus GenBank KU749412.1, with the POL protein of the Visna virus GenBank L06906.1
IKKKDSTKWRKLVDFRELNKRTQDFWEVQLGIPHPAGLKKKKSVTVLDVGDAYFSVPLDESFRKYTAFTI--PSTNNETPGIRYQYNVLPQGWKGSPAIFQSSMTKILEPFRAQNPEIVIYQYMDDLYVGSDLEIGQHRAKIEELREHLLRWGFTTPDKKHQKEPPF--------LWMGYELHPDKWTVQPIQLPE--KDSWTVNDIQKLV
|||| | ||| | ||||||| | |  | ||| ||| ||  || || || |||||  || |  | || |    |  ||  |  ||   ||||||| |||  |  | |||       | |    |||| | ||||  |      || |                                | |||| ||  |   |||      | |  ||||
IKKK-SGKWRMLIDFRELNKQTEDLAEAQLGLPHPGGLQRKKHVTILDIGDAYFTIPLYEPYRQYTCFTMLSP--NNLGPCVRYYWKVLPQGWKLSPAVYQFTMQKILRGWIEEHPMIQFGIYMDDIYIGSDL--GL-----EEHRGIVNELASYIAQYGFMLPEDKRQEGYPAKWL-GFELHPEKWKFQKHTLPEITEGPITLNKLQKLV

Are these matches,

By chance or convergence,
HIV1 and Visna virus share a recent common ancestor or,
the common ancestor was closer Visna virus from where HIV descended?

Could it be that HIV1 really comes from some kind of experiment involving Visna virus, as I read on the Internet?
The disease caused by the Visna virus is (100% deadly see https://pubmed.ncbi.nlm.nih.gov/16425963/ ). Could there be parallels to sars2?
But look at it: the Visna virus is transmitted by air the first symptoms are like pneumonia (or a brain infection) and otherwise are chronic progressing slowly like SARS-CoV-2 infection resulting in long Covid.
Symptoms are also felt after many days that the virus has come into action, this seems to me not to happen for coronaviruses or for flu viruses but for HIV and Visna virus.
Related links:

https://pubmed.ncbi.nlm.nih.gov/8289690/

https://pubmed.ncbi.nlm.nih.gov/16425963/

(preprint sars2 as retrovirus) https://www.biorxiv.org/content/10.1101/2020.12.12.422516v1

https://www.gao.gov/assets/100/91319.html

one 1966 HIV1 virus in https://www.pnas.org/content/117/22/12222  MN082768

They are both lentiviruses and share a distant common ancestor.
HIV-1 and HIV-2 are descendents from simian immunodeficiency virus (SIV). The following tree shows the relationships very clearly, from Wertheim and Worobey (2009) Dating the Age of the SIV Lineages That Gave Rise to HIV-1 and HIV-2 PLoS Comp Biol. here. The evolution of HIV-1 is heavily mixed with chimanpzee (cpz) SIV, and there are regions of highly anomalous homology juxtaposed with regions of notably lower homology.

and where it came from "SIVcol" and what is his genome?

Colobus monkey https://retrovirology.biomedcentral.com/articles/10.1186/1742-4690-10-107

so you say visna->siv->HIV1 & HIV2?  what about the chain is reversed visna->HIV1 & HIV2 ->SIV

No, he's saying OriginalLentivirus->SIV/HIV and OriginalLentivirus->visna.

So seems to me understand genome SIVcol can be find in GenBank with number
AF301156.1 (SIVcolCGU1), with other SIVcol gag protein KF214242.1,
with other SIVcol partial genome KF214241.1, the one base to all SIV Tree.
SIVcol the one all other SIV they claim came from, right?

@mgkrebbs no in the link they claim siv->HIV1 & HIV2  so first was siv virus than HIV1 and HIV2 they say

@RosLuP, yes that's true too, and well known, but the point is there is no evidence that HIV-1 is anymore related to visna than to other lentiviruses, such as the lentiviruses of cattle and the cat family.