cannin/biostars_qa huggingface dataset:
Question: I'm trying to decide what would be the most appropriate threshold criteria? In particular, this is for PFAMs and a human metagenomics dataset. I'm curious what people have done in the past? It looks like anvio uses --cut_ga as the default. https://github.com/merenlab/anvio/issues/498

**What been generally accepted as the most appropriate for metagenomics?**

--cut_ga
Use Pfam GA (gathering threshold) score cutoffs. Equivalent to
--globT --domT , but the GA1 and GA2 cutoffs are read
from the HMM file. hmmbuild puts these cutoffs there if the
alignment file was annotated in a Pfam-friendly alignment format
(extended SELEX or Stockholm format) and the optional GA annota-
tion line was present. If these cutoffs are not set in the HMM
file, --cut_ga doesn't work.

--cut_tc
Use Pfam TC (trusted cutoff) score cutoffs. Equivalent to
--globT --domT , but the TC1 and TC2 cutoffs are read
from the HMM file. hmmbuild puts these cutoffs there if the
alignment file was annotated in a Pfam-friendly alignment format
(extended SELEX or Stockholm format) and the optional TC annota-
tion line was present. If these cutoffs are not set in the HMM
file, --cut_tc doesn't work.

--cut_nc
Use Pfam NC (noise cutoff) score cutoffs. Equivalent to --globT
--domT , but the NC1 and NC2 cutoffs are read from
the HMM file. hmmbuild puts these cutoffs there if the alignment
file was annotated in a Pfam-friendly alignment format (extended
SELEX or Stockholm format) and the optional NC annotation line
was present. If these cutoffs are not set in the HMM file,
--cut_nc doesn't work.

http://www.cbs.dtu.dk/cgi-bin/nph-runsafe?man=hmmsearch

Answer: As explained above, none of these options will matter unless your HMM has the cutoffs set. The appropriate lines in HMMs look like this:

GA 25.00 25.00;
TC 25.00 25.00;
NC 24.90 24.90;

These bit-scores have to be set manually during model building, and their main purpose is to catch borderline matches that may be true hits but have statistically insignificant E-values. On rare occasions, the use of these cutoffs would eliminate a match that has barely significant E-value. A simple way of setting these scores is to pick the worst score in a known group of trusted family members.

Most of the time you will get the same result from these scores or from E-values, and I use the latter as a guide. It is good practice to manually inspect hits that are just above or just below the E-value threshold, and to set the database size (`-Z` switch) to a fixed number so that searches done over large period of time and with different databases can be compared. If you check [**how Pfam does it**][1], you will notice that they are not using cutoffs even though the HMMs support them. Instead, they set the database size to `45638612`.

[1]: https://pfam.xfam.org/family/PF00041#tabview=tab6
Question: I would like to create a LD vs distance(cm) plot in R using an output from PLINK. I have been trying with different --ld-window-kb numbers but I can not plot them in R. Furthermore, I have a question about average r2. In some papers ld_decay is plotted average r2 vs distance and some papers used r2 vs distance. Can everybody explain more?
thanks for your attention

Answer: For LD block length distribution, I ran this on the thinned SNP set.

plink --bfile "snp-thin" --blocks --blocks-max-kb 200 --out "snp-thin"

This produces snp-thin.blocks.det. Then some R.

dfr <- read.delim("snp-thin.blocks.det",sep="",header=T,check.names=F,stringsAsFactors=F)
colnames(dfr) <- tolower(colnames(dfr))

# ld block density
p <- ggplot(dfr,aes(x=kb))+
geom_density(size=0.5,colour="grey40")+
labs(x="LD block length (Kb)",y="Density")+
theme_bw()

ggsave("snp-thin-ld-blocks.png",p,height=8,width=8,units="cm",dpi=250)

![enter image description here][5]

Now you see LD block length distribution for your set of samples.

[5]: https://image.ibb.co/mjzPOH/snp_thin_ld_blocks.png
Question: I went to [the blast ftp database][1], there are 18 `nt` files, each is less than 800 MB, and for `refseq_genome` it has 83 files, most of which are larger than 800 MB, which means the `refseq_genome` is much larger than `nt` database. However, when I search the definition of `nt` on http://www.ncbi.nlm.nih.gov/BLAST/blastcgihelp.shtml, it says `nt` database include All GenBank + RefSeq Nucleotides + EMBL + DDBJ + PDB sequences (excluding HTGS0,1,2, EST, GSS, STS, PAT, WGS). No longer "non-redundant".

My question is:

1. In my understanding RefSeq Nucleotides should include `refseq_genome` and `refseq_rna`, so `refseq_genome` should be much smaller than `nt` database. why is `refseq_genome` alone is much larger than the whole `nt` database?
2. I tried one accession number `NZ_AARG01000001.1` from refseq bacteria genome, and blastn against `nt` and `refseq_genome` database. For `nt` case, it took a few seconds and got less than 10 hits. For `refseq_genome` database, it took more than 10 minutes and got more than 100 results (all the accession number began with NZ). Then I searched NZ and found NZ represent not completed project. So the difference between `nt` and `refseq_genome` is that nt doesn't include NZ records?

[1]: ftp://ftp.ncbi.nlm.nih.gov/blast/db/

Answer: On ["RefSeq accession numbers and molecule types"][1] you will see that the RefSeq accessions with the prefix `NZ_` are from whole genome shotgun (WGS) projects. As such these are excluded from 'nt'. Looking through the other 'genomic' sections of RefSeq, many of these are from WGS projects and are thus also excluded.

From the [NCBI BLAST][2] pages 'nt' is currently defined as:

> Title:Nucleotide collection (nt)
> Description:The nucleotide collection consists of GenBank+EMBL+DDBJ+PDB+RefSeq sequences, but excludes EST, STS, GSS, WGS, TSA, patent sequences as well as phase 0, 1, and 2 HTGS sequences. The database is partially non-redundant. In some cases identical sequences have been merged into one entry, while preserving the accession, GI, title and taxonomy information for each entry. Merged sequences include GenBank and RefSeq entries with identical sequences. Sequences added to the database since April, 2011 have also been merged with identical existing entries.
> Molecule Type:mixed DNA
> Update date:2014/07/24
> Number of sequences:23840180

In contrast 'refseq_genomic' is defined as:

> Title:NCBI Genomic Reference Sequences
>Molecule Type:mixed DNA
> Update date:2014/07/23
> Number of sequences:6733817

Note the difference in the number of sequences. However the 'refseq_genomic' is much larger when you look at the number of bases: 435,293,002,525 vs. 62,649,172,490. This is due to 'refseq_genomic' including assembled contigs, and whole chromosome assemblies, which are excluded from 'nt'.

[1]: http://www.ncbi.nlm.nih.gov/books/NBK21091/table/ch18.T.refseq_accession_numbers_and_mole/
[2]: http://blast.ncbi.nlm.nih.gov/
Question: This might be an XY question, so I'll explain my premise:

1. I have 3 VCF files, `f1`, `f2` and `f3`.
2. `f1` is an annotated VCF covering 50 samples
3. `f2` is an annotated VCF covering 5 samples, but only sites that are not in `f1`
4. `f3` is an un-annotated VCF covering 5 samples across sites in `f1` as well as not in `f1`
5. All annotations are site-level

I now wish to get this as one VCF files with all sites annotated and all sample-level information present.

When I merge `f1` and `f2`, I get a VCF with all annotated sites and all samples, but for those sites overlapping with `f3`, the `GT/AD/...` fields are empty, because that information is in `f3`. How do I merge these three datasets?

####Question:
In essence, can I do an operation to update genotype fields in one VCF file based on a sample+site match in another VCF file? If they were 2 `data.frame`s, the operation would be something like `vcf1[site, sample] <- vcf2[site, sample]`.

#### Current solution:

The way I see it, I might have to subset `f3` to `f1`-sites only, then `bcftools merge >` - that way I do not add any site, only samples. Then I `bcftools concat > `, so this time I add only sites, no samples. Any other solution will be appreciated. That solution does not work as `bcftools concat` cannot work on VCFs with different samples in them.

Answer: Here's my current solution:

1. Subset all `f1`-sites present in `f3`: `bcftools isec -n=2 -w1 -c none -o f3_subset f3 f1`
2. Pull annotations from `f1` into `f3_CommonSites_subset`: `bcftools annotate -c INFO -a f1 -o f3_subset_anno f3_subset`
3. Concat the new annotated file with `f2` to get all site annotations for the 5 samples: `bcftools concat -o f3_plus_f2 f2 f3_subset_anno`
4. Merge f1 and this 5-sample file to get final VCF: `bcftools merge -m none -o final_vcf f1 f3_plus_f2`

Just realized while I was writing this, I could just do `bcftools merge -Ou -m none f1 f3 | bcftools annotate -c FORMAT -a f2 -o final_vcf -`, so that way I would pick up the FORMAT fields exactly as I intended in the first place.

If anyone has a better solution, please add it in! Thank you!
Question: Hello,

I am trying to access information on genomic 3'UTR end and start positions using the Ensembl Biomart tool. However, there appear to be some transcripts which are labelled with the same transcript ID but nonetheless have different 3' UTR annotations.

e.g.:

```
Transcript ID 3' UTR start 3' UTR end
ENST00000474604 32793169 32793300
ENST00000474604 32792445 32792726
ENST00000474604 32791848 32791958
ENST00000474604 32791565 32791596
ENST00000474604 32790888 32791376
```

I was just wondering if there was anyway to get Biomart (or any other Ensembl tool) to print out the stable ID version increment (e.g. ENST00000474604.1, ENST00000474604.2 etc.)? Otherwise, trying to use this data it is going to be a bit of a hassle.

Also, I have downloaded the full cDNA set for Homo sapien (GRCh38) for processing and this includes the version increment - Ideally, I would like to easily map the data from the cDNA fasta file to the genomic co-ordinates obtained from BioMart

I have combed through biomart to see any option for including transcript version number but I can't seem to see anything. I am confused as to why this information would be omitted from the output

Thanks

EDIT: My question lead from an interpretation of Biomart output which was based on a misconception. See Sean Davis's post for more details

Answer: I think the understanding of your results is perhaps not quite right (or I am misunderstanding your question). The five regions that you give are not for different versions of the transcript. They signify the fact that this particular transcript has five 3'-UTR exons. If you want the UTR start and end, you can take the minimum and maximum of the two columns; which is which will depend on the strand of the transcript.
Question: Hi,

I'm looking for any technical information regarding how subread-align distinguishes properly paired from not-properly paired alignments for paired-end data (example summary below). I've looked through the subread/Rsubread documents and not found it (although I could have missed it).

Further, if anyone has experience aligning reads to a congener species assembly with subread and suggestions for parameter optimization, input would be greatly appreciated.

Thanks!

//================================ Summary =================================\\
|| ||
|| Total fragments : 12,563,751 ||
|| Mapped : 8,918,955 (71.0%) ||
|| Uniquely mapped : 8,918,955 ||
|| Multi-mapping : 0 ||
|| ||
|| Unmapped : 3,644,796 ||
|| ||
|| Properly paired : 5,290,912 ||
|| Not properly paired : 3,628,043 ||
|| Singleton : 2,281,674 ||
|| Chimeric : 80,494 ||
|| Unexpected strandness : 25,231 ||
|| Unexpected fragment length : 1,225,432 ||
|| Unexpected read order : 15,212 ||
|| ||
|| Indels : 145,958 ||
|| ||
|| Running time : 7.7 minutes ||
|| ||
\\============================================================================//

Answer: See arguments `minFragLength`, `maxFragLength` and `PE_orientation`. I think "Properly paired" means that both ends map to the same chromosome with the expected orientation and the fragment defined is within the expected length range.
Question: I have a file from Chip-Seq data. I want to use this information to look for the GC content of promoters.

Which software I can use?

Thanks
Ankur

Answer: The simplest would be to use nuc option in bedtools.

bedtools nuc -fi $ref-fasta-file -bed $promoterbedfile > $outputfile

Similar questions have been answered before -

- https://www.biostars.org/p/47047/
- https://www.biostars.org/p/70167/
Question: I am wondering what do the pieces of black lines between the coverage and reads mean? Thanks if anybody knows

![enter image description here][1]

[1]: /media/images/d3699f66-d82d-49f3-b90f-6475db16

Answer: Hi. They are this.

Found it on --> https://software.broadinstitute.org/software/igv/AlignmentData

> Downsampled reads areas are marked with a black rectangle just under the coverage track. The coverage track represents coverage for all the reads.

If you wanna get rid of the, go to preferences and turn off downsampling.
Question: Hi,

I have been trying to find the difference between the above two online for a while now, but I haven't got a satisfactory answer. I also didn't find a similar question on Biostars, so I thought of formally asking it now.

Tximport (and maybe other tools too) gives a couple of outputs for each gene, and two of them are - **abundance** and **counts**. What is the difference between them?

[This paper][1] gives a general idea that count based methods assign reads to genes directly, whereas abundance based methods assign abundance of each transcript with a probabilistic model that makes use of info such as fragment length distribution etc.

So, having said that, is this really the difference between the abundance and count values that I get for any gene from Tximport (or any tool in general)? And, in which situation is one of them a more meaningful/desirable quantity?

[1]: https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1004393

Answer: A count is simply that, a count of reads on some feature. An abundance is a more biologically meaningful (though not necessarily statistically useful) quantification of expression of a gene or transcript that is normalized in some way. Most commonly in this is TPM or some variant of that, but it could also be "copies per cell", which would be an abundance metric you could get from rt-qPCR. In other words, normalized counts aren't an abundance estimate since reads aren't a thing present in the cell, but an artifact of how we perform library prep and sequencing. The exception to this would be if you use a minion or equivalent to sequence full-length transcripts, since then a normalized count would estimate the abundance (on some likely relative scale) of a transcript in a cell or tissue.
Question: Dear community,
i have trees (<3000) in newick format with four species like this example:

((Spec4:0.529207,(Spec3:0.0803395,Spec2:0.0124315)),Spec1:0,Spec1:0);

I am only interested to detect the trees in which two species are clustering together, like in the example Spec3 and Spec2.
Is it possible to do that with a simple script or does anybody knows a software (tried phybin, ete3 compare already).
I will be grateful if you someone could help.

Answer: I think this might work, but it's a sort of 'brute force' way to do it. I would maybe re-factor your trees to cladograms and remove the branch lengths via a regex for the branch length and colon (in whatever your favourite regex language is), then you could simply `grep` or string search in some other manner for `(Spec3,Spec2)` and you'll find all trees which contain that grouping pretty easily.

e.g.: Remove decimals, sole zeros and colons from the file (probably not the most elegant regex):

Given your tree:

((Spec4:0.529207,(Spec3:0.0803395,Spec2:0.0124315)),Spec1:0,Spec1:0);

One could do:

cat test.tree | sed -e 's/[0-9]*\.[0-9]*//g' -e 's/0//g' -e 's/://g'

Yeilding:

((Spec4,(Spec3,Spec2)),Spec1,Spec1);

Then you can string search your yielded trees:

egrep -r -l "Spec(2|3),Spec(2|3)" .

Will give you all the filenames where Species 3 and Species 2 are adjacent nodes (in either orientation).

If you want to keep branch length in your trees as you're not just interested in topology, you could concoct a regex for use with `grep`:

egrep "Spec(2|3):(0?|[0-9]+\.[0-9]+),Spec(2|3):(0?|[0-9]+\.[0-9]+)" treefile.tree

But having to conjure that regex for every possible combination of topologies looks awful to me, so I'd be inclined to try it without the branch lengths.

I don't know how many topologies you're interested in finding in all your trees - this approach may not be feasible if it's a prohibitively large number.

-----

Slightly more complex, if you'd like to see the match, and the file name, this is an option:

2 example sed-treated trees:

((Spec4,(Spec5,Spec6)),Spec2,Spec3);
((Spec4,(Spec3,Spec2)),Spec1,Spec1);

Passing a 'dummy filename' in the form of `dev/null` tricks grep in to printing the filename (as it thinks it's working on multiple files) and the actual match itself by default:

for file in *.tree ; do egrep "Spec(2|3),Spec(2|3)" "$file" /dev/null ; done

Would yeild:

sed2.tree:((Spec4,(Spec5,Spec6)),Spec2,Spec3);
sed.tree:((Spec4,(Spec3,Spec2)),Spec1,Spec1);

With the appropriate string matches highlighted (if your terminal is configured for it).
Question: I'm doing a GWAS using ~15 million variants and ~800 people. I am unfamiliar with Linux, so I have tried using PLINK MDS and PCA functions to obtain principal components to be used as covariates in the association analysis to control for population stratification. When I plotted the p-values (QQ plot) obtained from the association analysis, the distribution was pretty messy, suggesting that I did not adequately control for population stratification. I took the following steps:

1. Pruned based on LD using PLINK --indep
2. Created a genome file:

./plink --bfile file --genome --extract plink.prune.in
3. Used --pca to generate an eigenvec file containing PCs

./plink --bfile gendep_merged --cluster --pca header --extract plink.prune.in --read-genome plink.genome

4. Performed the association analysis using 10 PCs from the eigenvec file as covariates:

./plink --bfile file --pheno phenotype.txt --allow-no-sex --covar plink.eigenvec --covar-name PC1,PC2,PC3,PC4,PC5,PC6,PC7,PC8,PC9,PC10 --out association --linear --adjust

Am I missing a step or should any of the flags used by modified in order to produce PCs that will adequately control for population stratification in this sample?

Any input would be greatly appreciated.

Answer: How exactly is using the first ten principle components controlling for "population stratification"? If I understand correctly, you're performing an association test, and telling the model fit to smooth out the ten biggest drivers of variance in your dataset? When you checked the principle components, did they indicate that the first ten explained the difference in population? Could you be smoothing out the effect you're testing for instead?
Question: I have a batch of 'noteworthy positions' (50,000+) and I'm trying to determine the best way to batch annotate them with my annotated genome (GTF format or anything which it can be converted to).

Does anyone know a package that already does this? I have considered writing a basic script but perhaps such a thing exists already.

My intention is to provide the coordinates and return columns such as:

```
Nearest Is within Intron/exon Distance positional annotation
gene gene/UTR/ from TSS (if exists)
Promotor?
```

Thanks

Answer: To get the nearest gene, you could use BEDOPS *closest-features*:

$ closest-features --closest positions.bed genes.bed > answer.bed

To answer if within gene/UTR/promoter/intron/exon, you could use BEDOPS *bedmap --indicator* (or similar options, depending on what level of detail you need):

```
$ bedmap --echo --indicator positions.bed genes.bed > answer_genes.bed
$ bedmap --echo --indicator positions.bed UTR.bed > answer_UTR.bed
$ bedmap --echo --indicator positions.bed promoter.bed > answer_promoter.bed
```

etc.

To get the distance from the nearest TSS, you could use BEDOPS *closest-features* again:

$ closest-features --closest --dist positions.bed TSS.bed > answer.bed

To get annotations of any overlapping elements, you could use BEDOPS *bedmap --echo-map*:

$ bedmap --echo --echo-map positions.bed annotations.bed > answer.bed

If your annotations are in a non-BED format, you could use BEDOPS *convert2bed* to turn them into BED:

```
$ gtf2bed < annotations.gtf > annotations.bed
$ gff2bed < annotations.gff > annotations.bed
$ vcf2bed < annotations.vcf > annotations.bed
```

etc.

Because BEDOPS tools use Unix I/O streams, you can convert inside a larger BEDOPS operation pipeline:

$ gff2bed < annotations.gff | closest-features --closest --dist positions.bed - > answer.bed

Or you can use process substitution, if you use `bash`:

$ closest-features --closest --dist positions.bed <(gff2bed < annotations.gff) > answer.bed

To turn BED-formatted annotations into BED-formatted lists of UTRs, introns, exons and promoters, etc. is easy and a cursory search of Biostars or seqanswers can help locate small scripts in *awk* or Perl etc. that can do this for various annotation subcategories.

For instance, to turn a BED-formatted list of stranded gene annotations into 1k proximal promoter regions:

```
$ awk '{ \
if ($6=="+") { \
print $1"\t"($2 - 1000)"\t"$2"\t"$4"\t"$5"\t"$6; \
} \
else { \
print $1"\t"$3"\t"($3 + 1000)"\t"$4"\t"$5"\t"$6; \
} \
}' genes.bed \
> promoters.bed
```

Then it is a matter of using `promoters.bed` in downstream BEDOPS operations.
Question: Hi everybody!

I submitted a list of gene to DAVID Functional annotation chart tool. I decided to used EASE (p-value) cutoff of 0,05 (from the DAVID publication on Nature Protocols in 2009) and selected Benjamini correction. Now I have problems with results interpretation. Is the Benjamini value the (i/m)Q value? So if it is bigger than my p-value, is my term significantly enriched? Or should I interpret the Benjamini value like a new p-value corrected already and is it significant if it is <0.05? (It looks strange for a multiple comparison tests).

Thank u

Answer: Benjamini correction is your false discovery rate and it is your adjusted p-value. So you should forget about your p-value after correction. So your test is significant if your adjusted p-value is smaller than criteria (such as 0.05 or 0.01). If you want to more about multiple testing, you can check [here][1]

[1]: https://en.wikipedia.org/wiki/Multiple_comparisons_problem
Question: Hi All,

I just started working with methylation (WGBS) data. I have used Bismark and it generated the methylation of Cs and consecutive Gs. I am guessing this Gs methylation is from Cs from other strand (please correct me if I am wrong). My question is should I consider both while reporting the methylation in a region or should I filter only Cs rows and perform the further analysis?

Thanks.

Answer: As you surmised, the "G" is actually the C on the - strand. The methylation of a region should include both strands, unless you know that whatever you're studying is only affected by single-stranded methylation. Note, however that the Cs in a CpG (the C and subsequent G in bismark's files) typically show roughly symmetric/equal methylation ratios, so a common strategy is to simply combine them into CpG-level methylation levels. Bismark likely has a method for doing that (if not, you can do it with MethylDackel).