cannin/biostars_qa huggingface dataset:
Question: Hi
This is probably a question for stack overflow but I am posting it here because there I am not getting much replies, apologies for this.

I have a 24 data frames with different number of rows and with various columns (for instance column SNPs ID and other columns with corresponding values for each SNPs) and I stored those data frames inside a list.
I am doing various operations on the data frames at the same time.
For instance if I want to order in decreasing manner a column in all the data frames inside the list I do:

myfiles_ordered<-lapply(myfiles, function(x) { x[ order(x$column_name_to_order, decreasing=T),]})

Now, after ordering that column I would like to take the top 5 % of the values of it.
I was thinking I can subset all the data frames based on their specific row number multiplied by 0.05 and I wrote something like this:

myfiles_top5<-lapply(myfiles_ordered, function(x) {x[1:nrow(x)*5/100,]})

However, it does not seems to work.
Any help highly appreciated, thanks.

Answer: Try:

myfiles_top5 <- lapply(myfiles_ordered, function(x) { x[ 1:round(nrow(x)*5/100), ]})

Because we are creating a sequence of 1 to n, then applying 5% for all of them.

1:nrow(mtcars) * 5/100
# [1] 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80
# [17] 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60

Not what we need...

Instead we need to get 5% then get the sequence, using parenthesis `()`:

1:(nrow(mtcars) * 5/100)
# [1] 1

Again not ideal as, below both give 1:

1:1.2
# [1] 1
1:1.6
# [1] 1

Whereas we might need `1:2` for `1:1.6`, so we use `round`:

1:round(1.2)
# [1] 1
1:round(1.6)
# [1] 1 2

--------

**Update:** We can do ordering and subsetting 5% in one go, e.g.:

# Using base
head(mtcars[ order(mtcars$mpg, decreasing = TRUE), ], round(nrow(mtcars) * 5/100))
# mpg cyl disp hp drat wt qsec vs am gear carb
# Toyota Corolla 33.9 4 71.1 65 4.22 1.835 19.90 1 1 4 1
# Fiat 128 32.4 4 78.7 66 4.08 2.200 19.47 1 1 4 1

# Using dplyr
library(dplyr)
top_n(mtcars, round(nrow(mtcars) * 5/100), wt = mpg)
# mpg cyl disp hp drat wt qsec vs am gear carb
# 1 32.4 4 78.7 66 4.08 2.200 19.47 1 1 4 1
# 2 33.9 4 71.1 65 4.22 1.835 19.90 1 1 4 1
Question: I have found out the core ortholog set (Core Genome) of 30 bacterial strains using NCBI Blast Package.
But finding Pan Genome (Unique genes + accessory genes + Core genes) of same dataset of 30 organisms is becoming hectic.
As it is not possible to align each genome with other, around overall 900 times.
I can`t derive any other logical pattern so that I can determine the accessory genes without repeats.
Please help.
Thanks in advance.

Answer: You probably would not need to align every gene against every other (if that's what you meant above). A simpler technique would be to add keep adding genes to a database if they are sufficiently different from the genes that are already there. As you perform the alignments you will need a simple program to tabulate which genes have been hit just once over the entire process (these are the unique), genes that were hit for every strain (core genes) the rest are the accessories.

You might also want to consult what the literature says:

The Pangenome Structure of Escherichia coli: Comparative Genomic Analysis of E. coli Commensal and Pathogenic Isolates
Question: How can I extract mutations/variants from Pubmed/PMC articles using terms and without need to download these articles?
Any command line tool I can do that?

Answer: NCBI Elink from a pubmed id (e.g: 21378989 ) maps to a dbsnp id (eg. rs387906747 ) .

wget -q -O - "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/elink.fcgi?dbfrom=pubmed&db=snp&id=21378989"

pubmed

21378989

snp
pubmed_snp

387906747
Question: Hey,

is there a way to tell which of the transcripts for a given gene corresponds to the canonical isoform in the Gencode dataset?
I couldn't find such information anywhere.

Best

Answer: [APPRIS][1] have carried out principal isoform analysis on the GENOCDE transcripts. These are also available in [Ensembl][2].

[1]: http://appris.bioinfo.cnio.es/#/
[2]: http://www.ensembl.org/index.html
Question: When using the [SPAdes Genome Assembler][1], you can use either the BayesHammer read error correction tool for Illumina reads or the IonHammer read error correction tool for IonTorrent data.

I want to use SPAdes in a general workflow and I want to know if there is an automatic way to determine if reads are either from Illumina or IonTorrent (without asking the user) so that I can select the correct read error correction tool.

[1]: http://bioinf.spbau.ru/en/spades

Answer: To determine if the reads are from illumina or IonTorrent, I think the best solution is to look at the read names format in the fastq files.

**Illumina reads :**

```
@M02681:36:000000000-AAUMD:1:1101:17649:1746 1:N:0:1
@:::::: :::
```

**Ion Torrent reads **:

@I7HE0:00523:00720
Question: I have been struggling with this for a while. I appreciate if someone can help:

I am running a linear model using limma package on a set of micro-array probes and a disease phenotype adjusting for a confounding factor. Ideally I want my model to be like this:

Disease phenotype (continuous) ~ Micro-array probes + confounding factor

But in Limma package, There is no way to have the probe matrix in the middle. It's always the response variable. Whatever I put in my model matrix will act as the independent variable as it is seen here, which is not what I want:

lmFit(probe_matrix, design=model.matrix(~Disease phenotype+confounding factor))

Is there a way around this?

Answer: I think I figured it out myself. Apparently this is the only way:

pvals<-sapply(1: dim(probe_matrix[1], function(x) eBayes(lmFit(Disease_phenotype, design=model.matrix(~t(probe_matrix[x,])+t(confounding factor))))$p.value)
Question: How can I filter a vcf filter a VCF file on minimum genotype depth and genotype quality for each sample.

I am looking for a way to filter variants from a VCF file by checking that all samples for a site pass 2 critera

sample.DP > 10
sample.GQ > 15

I thought I could do this using bcftools

bcftools view -i 'MIN(FMT/DP>10) && MIN(FMT/GQ>15)' my.vcf.gz

Somehow this include expression does not seem to be applied.

In my output are still lots of variants with genotypes likes

GT:AO:DP:GQ:PL:QA:QR:RO
0|0:0:16:3:0,48,502:0:554:16 (=depth 16, quality 3)
1|1:1:1:3:21,3,0:37:0:0 (=depth 1, quality 3)

I tried putting the include commands in two different bcftools commands

bcftools view -i 'MIN(FMT/DP>10)' | bcftools view -i 'MIN(FMT/GQ>15)' my.vcf.gz

And I tried without the FMT prefix for the genotype quality

bcftools view -i 'MIN(FMT/DP>10) && MIN(GQ>15)' my.vcf.gz

Still I get variants back with genotypes that don't match the criteria.

Am I misunderstanding the bcftools include expression documentation?

https://samtools.github.io/bcftools/bcftools.html#expressions

Or is there another way that I can achieve this filtering step?

Answer: The problem was wrong place of the parentheses:

Does not work as expected:

bcftools view -i 'MIN(FMT/DP>10) & MIN(FMT/GQ>15)' my.vcf.gz

Works as expected:

bcftools view -i 'MIN(FMT/DP)>10 & MIN(FMT/GQ)>15' my.vcf.gz
Question: Hi all,

I checked the Nanopore basecalled reads using MinION QC tool. From what I've understand, I need to select only the reads which are passed the N50 value. MinION QC gives the information only. Now how to select those passed reads from the data. Which tool to be used to select the best quality reads and to trim?
I am not a bioinformatician. But trying to become one.

Answer: You may be [able to use][1] `nanofilt` which is part of a suite of tools called `nanopack` written by @Wouter DeCoster. If you have not used these tools before you may want to get the [full suite][2].

Note: Trimming and filtering are two different operations. Since title of your post says trimming, you may be able to use other trimmers (`bbduk.sh` from BBMap suite, trimmomatic etc) for that analysis.

[1]: https://github.com/wdecoster/nanofilt
[2]: https://github.com/wdecoster/nanopack
Question: Hi all,

I downloaded a RNAseq raw data and found it in weird format. The data is generated by ABI_SOLID platform. Data ID is SRR179593.fastq and here are a few lines inside it:

@SRR179593.1 mendel_20110320_FRAG_BC_Ryan_RNA_Seq_2_58_240/1
T.20.2312.0100.0111200312.0.1300.0.01.0202.0.100020
+
!!<:!7<=)!>/?:!78;-.A6-3%!&!(*?%!0!%,!.*%9!0!)5;1%+
@SRR179593.2 mendel_20110320_FRAG_BC_Ryan_RNA_Seq_2_58_300/1
T.30.2010.0313.1231232001.1.1112.0.12.0210.0.120003
+
!!9%!>48'!9:.%!7;%058A2%:!5!%<>;!+!-.!+'*'!)!',,)&+
@SRR179593.3 mendel_20110320_FRAG_BC_Ryan_RNA_Seq_2_58_638/1
T.23.1101.3222.2301222112.2.1332.2.22.2212.1.002023
+
!!-'!'%3(!'%3'!1%).).1)%'!5!%2-7!'!'5!%%0-!%!(.'(%+

I tried to map it to the genome using STAR:

#!/bin/bash
cd .
/usr/local/star/2.4.1c/bin/STAR --runThreadN 4 \
--genomeDir /escratch4/tw30282/tw30282_Mar_17/STAR/ENSEMBL.homo_sapiens.release-75 \
--readFilesIn SRR179593.fastq.gz \
--outFileNamePrefix Wei_NPC_rep1_STAR_SE. \
--outSAMtype BAM SortedByCoordinate \
--readFilesCommand zcat

But an error returned and said input file in wrong format.

Is there anyway to transform or fix the format?

Thanks,
Tianming

Answer: ABI SOLiD data is in colorspace rather than basespace, which is why you have 0123 instead of ACGT for the sequence. You need to use an aligner that works with SOLiD data. Normally, if you weren't getting the data from SRA the SOLiD data would have 2 separate files, one with the sequence reads and one with the read qualities.

Tophat2 works with colorspace data, you have to specify the options that your data and genome index are in colorspace, but I don't see anything in the Hisat manual so I don't think it does.
Question: Hi everyone,

Is there a smart way to get all contig boundaries after assembly?
I know samtools can get number of sequences covering a defined region. Then I could write a pipeline to find each contig's boundary by looking for a region flank by coverage of 0. It seems this costs a lot of time.
Is there an alternative and quicker way to know the start and end position of each contig?

Thanks a lot.

Dadi

Answer: You could use samtools mpileup and awk to get the contigs:

samtools mpileup file.bam |\
cut -d ' ' -f 1,2 |\
awk -F ' ' 'BEGIN {chr="";start=-1;end=-1} {if(chr!=$1 || int($2)!=end+1) { if(chr!="") {printf("%s:%d-%d\n",chr,start,end);} chr=$1;start=int($2);end=int($2);} else { end=end+1;}} END {if(chr!="") {printf("%s:%d-%d\n",chr,start,end); } }'
chr1:10044-10300
chr1:10330-10468
chr1:10504-10557
chr1:10585-10633
chr1:11302-11355
chr1:11594-11647
chr1:11696-11823
chr1:11827-13107
chr1:13122-13840
chr1:13874-13973
(...)
Question: Hi! I have a PGS file with the weights for each of the variants for a specific disease, CAD. The goal is to now calculate PRS scores for individuals in the UK Biobank genetic data. I have the plink bed, bim and fam files from the UK biobank data. What would be the steps to prepare the "Target data" from the uk biobank data?

I understand that using the PGS file as the "base data" requires adding a "fake" p-value column. I have tried using the bed, bim and fam files from uk biobank along with a GWAS file as the base data to which it gives: "Error: All sample has invalid phenotypes!", "Errorr: No sample left" . What am I missing here?

> PRSice 2.3.5 (2021-04-06)
https://github.com/choishingwan/PRSice
(C) 2016-2020 Shing Wan (Sam) Choi and Paul F. O'Reilly
GNU General Public License v3
If you use PRSice in any published work, please cite:
Choi SW, O'Reilly PF.
PRSice-2: Polygenic Risk Score Software for Biobank-Scale Data.
GigaScience 8, no. 7 (July 1, 2019)
2021-06-02 18:41:16
./bin/PRSice \
--a1 a1 \
--a2 a2 \
--bar-levels 0.001,0.05,0.1,0.2,0.3,0.4,0.5,1 \
--base CAD_UKBIOBANK.gz \
--beta \
--binary-target F \
--bp bp \
--chr chr \
--extract PRSice.valid \
--interval 5e-05 \
--lower 5e-08 \
--no-clump \
--num-auto 22 \
--out PRSice \
--pvalue pval \
--seed 1999345467 \
--snp oldID \
--stat beta \
--target chr# \
--thread 4 \
--upper 0.5

Initializing Genotype file: chr# (bed)

Start processing CAD_UKBIOBANK
==================================================

SNP extraction/exclusion list contains 5 columns, will
assume first column contains the SNP ID

Base file: /shared/Jenish/CAD_UKBIOBANK.gz
GZ file detected. Header of file is:
uniqid chr bp a1 a2 beta se pval N af oldID info zval

Reading 100.00%
7947837 variant(s) observed in base file, with:
1202749 variant(s) excluded based on user input
6745088 total variant(s) included from base file

Loading Genotype info from target
==================================================

488377 people (223459 male(s), 264780 female(s)) observed
488377 founder(s) included

181798 variant(s) not found in previous data
602458 variant(s) included

There are a total of 1 phenotype to process

Processing the 1 th phenotype

Error: All sample has invalid phenotypes!
Error: No sample left

Error:
Execution halted

enter code here

Answer: You also need the `--no-clump --no-regress --bar-levels 1 --fastscore` parameters to tell PRSice not to perform clumping (as I supposed you are using pre-computed effect sizes, which has already accounted for LD) and not to perform the regression as you don't need to optimize the p-value thresholds

(I assume you are using PRSice as I recognize the error message)
Question: Hi,

I know that expression data generated by different experiments cannot be combined directly for downstream analysis due to "batch-effects". My question is very straight forward. Is batch effect removal necessary before merging normalized microarray datasets that use the same Affymetrix platform but different chip types, say HG-U133A and HG-U133PLUS2. Or simple merging can solve the purpose?

All kind suggestions will be appreciated.

Thanks

Answer: boxplot can easily figure out it.

you can use combat (http://www.bu.edu/jlab/wp-assets/ComBat/Abstract.html)
Question: I have created my own genome for research. I have data for this fake genome. Now I want compare this genome with a BAM file on UCSC or any other browser. However, the browser only allows me to select a known genome as the reference. How do I provide my own genome on a browser, it doesn't have to be the UCSC browser.

Answer: Hello,

[Galaxy][1] has a very robust "Custom Reference Genome" feature that functions with both analysis tools and visualization. This permits nearly any fasta file to be used as a reference genome. **Galaxy's visualization browser is called "Trackster".**

Start by loading the fasta file of your "fake" genome. From there it can be directly used with tools. To go forward and use this genome with Trackster, create a "Custom Build" based on the fasta version of the genome. All of this takes a few minutes to do at most. Full instructions, including help videos, are in the http://galaxyproject.org wiki here:
Summary: http://wiki.galaxyproject.org/Support#Custom_reference_genome
Details & Troubleshooting: http://wiki.galaxyproject.org/Learn/CustomGenomes

Hopefully this will work out for you, Jen, Galaxy team

[1]: http://usegalaxy.org
Question: I had my whole genome sequenced. In plink, I put the VCF file through and it produced 3 files: plink.bed, plink.bim, and plink.fam. As recommended by Sam in another comment, I made a backup of plink.bim and ran this Linux command on the backup: "sed 's/^chrM\s/25\t/g; s/^chrX\s/23\t/g; s/^chrY\s/24\t/g; s/^chr//g' plink-backup.bim > plink.bim" With this the changes are reflected in the file plink.bim. I then run admixture "./admixture -j2 plink.bed 4" and I got this error:

**

Error: detected that all genotypes are missing for a SNP locus.
Please apply quality-control filters to remove such loci.

**

I have no idea what to do at this point.

Answer: So it sounds like one of your SNPs has all alleles missing, which is likely to happen if you've got only one genome there (yours).

You can remove SNPs that have missing alleles using plink

plink --geno 0.1 --file your_file --recode --out your_filtered_file

geno 0.1 means 'remove variants with missing call rates exceeding the provided value (default 0.1)', 0.1 (10%) should be more than enough for your case since your missing rate is 1.0 (100%)

This should print how many SNPs have been filtered out.
Question: can someone explain for me that why there are a few gaps in reference genome??

thanks!

Answer: The simple answer is that certain stretches of a genome contain sequence that is difficult to sequence, mainly due to repetitive regions, tracks of the same base, GC composition, closed DNA, etc. Searching Google for "gaps genome" brings up a whole host of references about the causes of gaps and the attempts to close them.
Question: So I have a pretty simple question that I actually can't seem to be able to figure out. So far I have been using bedtools to find overlaps between two bed files using intersectBed. Example:

```
Bed A
Bed B
Bed C
```

Now I have three generated bed files that I want to overlap to find peaks/regions common among all three.

ABC

The multiintersect option doesn't seem to have any documentation besides the information found here https://www.biostars.org/p/13516/ but the function itself doesn't seem to give me the information I'm looking for. Specifically I want to feed in three bed files and find only the common regions between ALL three, not AB, AC, BC, and ABC in one large file which seems to be the output shown in the linked example.

I believe that `intersectBed -a A -b B C` does something similar to the above, but perhaps I'm simply running it wrong and my attempts are errors on my part.

Can this be done using bedtools? If not, what other similar software is out there that can accomplish this?

I appreciate any help,

Carlos

Answer: Answer:

First find in-between any two files and use the results to compare with third file.

intersectBed -a 1.bed -b 2.bed | intersectBed -a - -b 3.bed
Question: Hello,

We have RNAseq data on a sample from a few years ago that looked odd (possibly mislabeled), so we re-sequenced a bunch of candidates that could possibly have been that sample.

Is there any way to find out which of the re-sequenced samples the original sample is closest to? I could run PCA but these re-sequenced candidates and the original would show some sort of batch effect and that might interfere with how the samples stratify on a PCA plot. Is there any way I could account for that?

I'd appreciate any pointers. Thank you!

Answer: You can try the plotPCA function from DESeq2 after normalizing data with `vst` which has a mode `blind` that can be set to TRUE and FALSE. In FALSE, it will respect the design formula so the batch effect would be respected. Pearson as suggested might make sense in combination with hierarchical clustering to visualize the results (see `corrplot` package). Another PCA-like metric is the Multidimensional Scaling (`plotMDS` function) from `edgeR`. In any case the individual approaches should at least in part confirm each other before you make a decision.

Apart from those dimensionality-reduction techniques, it **might** be possible to perform variant calling (even though difficult from RNA-seq) and see if some characteristic SNPs or patterns of SNPs can give valuable information. You could take variants that are common in humans (e.g. find some with high MAF in the hg38 VCF from dbSNP latest release containing 1KG and TOPMED MAFs) and then see if for those variants the resequenced sample is similar to the odd ones.
Question: I have a case-control dataset and I want to perform logistic regression and conditional logistic regression based on HLA multi-allelic data, using r. I want to find the effect on specific alleles on a trait. How do I do this what is the format. Most examples are based on SNP biallelic data.
For instance at HLA-A I may have up to 30 unique alleles, at HLA-B it could be 50. Should I recode all the alleles and perform logistic regression on genotype pairs?

Answer: Find a way to produce a data frame containing the counts of each alleles that you see, and the case-controls status. E.g. (fake data)

> df

DX DRB1.0401 DRB1.0404 DRB1.0405 DRB1.0408
1 0 0 0 1 1
2 0 0 0 0 2
3 0 0 0 0 2
4 1 1 0 0 1
5 1 1 0 0 1
6 1 0 1 0 1

If you are interested on the effect of a specific allele, then you can do, e.g.

summary(glm( DX ~ DRB1.0401 , family="binomial", data=df ) )

If you are interested in the effect of your HLA locus as a whole, then you can do, e.g.

full<- glm( DX ~ DRB1.0401 + DRB1.0404 + DRB1.0405 + DRB1.0408,
family="binomial", data=df )
null<- glm( DX ~ 1 , family="binomial", data=df )

anova( null, full , test="Chisq")

adding covariates to the models if deemed necessary.
Question: If I want to read a bam like this

#include
#include
#include
#include

int main(int argc, char *argv[]){

samFile *fp_in = hts_open(argv[1],"r"); //open bam file
bam_hdr_t *bamHdr = sam_hdr_read(fp_in); //read header
bam1_t *aln = bam_init1(); //initialize an alignment

char *chrom = argv[2];
int locus = atoi(argv[3]);
int comp ;

printf("%s\t%d\n", chrom, locus);

//header parse
//uint32_t *tar = bamHdr->text ;
//uint32_t *tarlen = bamHdr->target_len ;

//printf("%d\n",tar);

while(sam_read1(fp_in,bamHdr,aln) > 0){

int32_t pos = aln->core.pos +1; //left most position of alignment in zero based coordianate (+1)
char *chr = bamHdr->target_name[aln->core.tid] ; //contig name (chromosome)
uint32_t len = aln->core.l_qseq; //length of the read.

uint8_t *q = bam_get_seq(aln); //quality string
uint32_t q2 = aln->core.qual ; //mapping quality

char *qseq = (char *)malloc(len);

for(int i=0; i< len ; i++){
qseq[i] = seq_nt16_str[bam_seqi(q,i)]; //gets nucleotide id and converts them into IUPAC id.
}

//printf("%s\t%d\t%d\t%s\t%s\t%d\n",chr,pos,len,qseq,q,q2);

if(strcmp(chrom, chr) == 0){

if(locus > pos+len){
printf("%s\t%d\t%d\t%s\t%s\t%d\n",chr,pos,len,qseq,q,q2);
}
}
}

bam_destroy1(aln);
sam_close(fp_in);

return 0;
}

How do I get the strand info? In sam files it seems to be represented as 16 for reverse, 0 for positive.

Answer: in sam.h:

/*! @abstract the read is mapped to the reverse strand */
#define BAM_FREVERSE 16

usage:

(aln->core.flag & BAM_FREVERSE)
Question: I have a very basic question. I have created a local (mysql) database in my system. I would like to know what would be the best option for a front end. I had used html before and want to try out something new. I might also add some perl code in between for some computational work. Kindly help me in this regard..

Answer: To let others to access and update the database, I recommend OpenOffice Base

- It is Open Source and cross platform
- It produces a nice(ish) user interface QUICKLY (if you have to make your own webpages it is a nightmare)
- On top of the server side permissions (who can UPDATE, DROP, SELECT...) you can put some other client side restrictions
- It allows to make, save and share custom QUERIES, FORMS (either to read or to update tables)

You can start here to see how to connect Base to mysql

I have just started using it. I am in a large(ish) group of biologist that will need to update different tables (about samples, DNA preparation and so on) so that such information can be used by me and another bioinformatician for analysis. Also people in the lab need to retrieve information and I think this is the way to go. Happy to help more (and share link on documentation that is sparse) if you go this way

If you want a nicer user interface (don't like mysql client) you can use MySQL Workbench

I use it mainly to keep an up to date EER model, but it is much more powerful than that