cannin/biostars_qa huggingface dataset:
Question: I'm trying to do imputation by IMPUTE2 with a VCF file.

To generate the input files, I used plink:
./plink --vcf file.vcf --recode oxford --chr 1 --out output

Then I got four files: output.gen, output.sample, output.nosex and output.log

The basic IMPUTE2 requires inputs including output.gen, and also a strand file indicating the strand information of each SNPs in the output.gen file.

I'm not very familiar with the chip data and I'm not sure if I understood correctly, but I remember all SNPs in VCF files should be on the + strand. Does that mean all SNPs in the .gen file converted from vcf are strand +?

How can I get this strand file for my data?

Thanks!

Answer: I found a tool to do this: https://github.com/biocore-ntnu/snpflip
Question: I'm using VEP to annotate variants in VCF files. In the analysis of the results I want to select only the protein coding transcripts. Should I look at the field "consequence" added by VEP?
I don't have a biological background and I didn't find any answer when searching on the internet until now.

Answer: I usually run VEP with option "--everything". In the results you can see the consequence for each overlapping transcript. You have to carefully parse the results, which is a bit complicate when there are multiallelic variants. There is a field for transcript-biotype named "BIOTYPE" and another field for the transcript ("Feature"). Just set a filter for BIOTYPE to be "protein_coding".

Alternatively, you may preload a list of protein-coding transcripts (you can get them from Biomart), and see whether the transcript in the "Feature" field is within your list.
BTW, I would also recommend not to use all the protein-coding transcripts but only the ones that are more reliable (e.g. with CCDS or APPRIS support).

[edited to correct a detail]
Question: Hi,

I have downloaded the current version of the stand-alone-blast (ncbi-blast-2.2.29+) and I am trying to use blast (blastn) to find similarity of of a group of nucleotide sequences that I have. However, I am interested on only the top 3 matches. I tried searching online and I saw some posts that suggests using `-K`, but I realized this does not work with the new version that I am using. I looked at the help document and I tried using (`-max_target_seqs`) and (`-num_alignments`) but none of them worked. The result contains all the matches found by blast.

Does anyone know how to limit the results to let say just top 3 matches?

Thanks!

Answer: Depends what you are trying to do.

As [Neilfws][1] says, if you want to limit the number of hits reported you can use (from the NCBI BLAST+ help output):

```
-num_descriptions =0>
Number of database sequences to show one-line descriptions for
Not applicable for outfmt > 4
Default = `500'
* Incompatible with: max_target_seqs
-num_alignments =0>
Number of database sequences to show alignments for
Default = `250'
* Incompatible with: max_target_seqs
```

These correspond to the '-v' and '-b' options in legacy NCBI BLAST:

```
-v Number of database sequences to show one-line descriptions for (V) [Integer]
default = 500
-b Number of database sequence to show alignments for (B) [Integer]
default = 250
```

The '-K' option in legacy NCBI BLAST:

```
-K Number of best hits from a region to keep. Off by default.
If used a value of 100 is recommended. Very high values of -v or -b is also suggested [Integer]
```

Is slightly different and maps to the '-culling_limit' parameter in NCBI BLAST+:

```
-culling_limit =0>
If the query range of a hit is enveloped by that of at least this many
higher-scoring hits, delete the hit
* Incompatible with: best_hit_overhang, best_hit_score_edge
```

You may also want to limit the number of matches reported per hit (i.e. limit the number of HSPs):

```
-max_hsps =0>
Set maximum number of HSPs per subject sequence to save (0 means no limit)
Default = `0'
```

For more information about the NCBI BLAST+ command-line options see:

- NCBI BLAST+ usage/help output (run the relevant BLAST program with the '-help' option)
- [BLAST Command Line Applications User Manual][2]
- [NCBI's BLAST help][3]

[1]: /u/66/
[2]: http://www.ncbi.nlm.nih.gov/books/NBK1763/
[3]: http://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs
Question: I want to know how a skipped region in the reference, or N in the CIGAR string, looks in the alignment. To try and explain what I mean I use the example provided from the SAM format specification (http://genome.sph.umich.edu/wiki/SAM), which does not include an N example:

```
Ref + read
RefPos: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
Reference: C C A T A C T G A A C T G A C T A A C
Read: ACTAGAATGGCT

Alignment
RefPos: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
Reference: C C A T A C T G A A C T G A C T A A C
Read: A C T A G A A T G G C T

Cigar:
POS: 5
CIGAR: 3M1I3M1D5M
```

---

Now, in position 11 there is an insertion in the reference sequence. However, I would think that you can't distinguish between a skipped region or an insertion in the reference. Therefore the CIGAR string could also have been `3M1I3M1N5M`

So how is it the alignment of a skipped region or an insertion in the reference sequence different? Is it only a skipped region if the C in position 11 is an N?

Answer: There's no *a priori* way to always distinguish between deletions (D CIGAR operations) and splicing (N CIGAR operations). In practice, most RNAseq aligners (e.g., tophat2 and STAR) have parameters with semi-arbitrary thresholds for the minimum intron size or maximum deletion size. In the case of STAR, any gap less than `alignIntronMin` (21 bases last I looked) is considered a deletion. I can't recall exactly what the tophat2 option for this is off-hand, but it's in there somewhere.

It's probably worth pointing out that the default values for these might be worth changing in some cases. I suspect that if someone were interested in splicing changes in cancer cells where there are a bunch of deletions that these parameters might need some tweaking (though presumably one would do WGS or WES alongside).
Question: Hi all,

I am working on microarray expression. I am new to this, so pardon me if I provide incomplete details.

I have downloaded the raw data from GEO (gilent-014850 Whole Human Genome Microarray 4x44K G4112F) and we want to look at expression levels of genes with different datasets added on top (expression levels of genes with and without epigenetic marks, transcription factors).

When I was extracting the information, I realised that there are multiple probes for certain genes. At first, I took the highest expression values for each gene and then later compared it with the average and the difference was notable.

I was wondering if there is a method to identify which probes are present in all isoforms and which are present in a few.

Thanks for the help.

Answer: Yes, I agree with mapping the probe sequences on your own. However, I typically wouldn't do this prior to analysis - in practice, there are a lot of probes to check.

I would typically conduct differential expression with all the probes, and allow a single probe to be sufficient for differenital expression (at which point, I could take some time to understand any significant differences between multiple probes that map to the same gene - for specific results that look interesting).

Your strategy should depend upon your method of integration. If you are looking for overlapping gene lists, I would recommend the strategy listed above. However, if you need to work with the absolute expression values in the integration process and you need a single expression value per gene, I would typically just average the expression among probes that map to the same gene.
Question: Hi,

Simple question but I'm a little bit confused due to different description I read.
So what is the difference between insert size and fragment size in DNA-Seq ( or RNA-Seq ) ?

Thanks

Answer: See the above figure ( from https://www.biostars.org/p/106291/ )
Question: Hello all,

I would like to know how to filter a VCF by number of individuals with genotype data (i.e. without missing data).

First, some background:

I have a large VCF file containing SNP data for 95 individuals called from bcftools mpileup output. BAMs passed to mpileup contained data for alignments of RNAseq reads to a de novo transcriptome. Thus, the VCF lists SNPs within each contig defined in the reference de novo transcriptome.

The 95 individuals were collected from 6 different populations. For population genetic analysis, I would like to produce a VCF that only contains SNPs at which there are data for at least 75% of individuals in each population (11 or 12).

I have a decent familiarity with bcftools, and know how to filter missing data, but I am not sure how to approach this. As expected, filtering on the N_SAMPLES variable alone doesn't seem to modify the VCF. Perhaps my best course of action would be to first split the VCF into population-specific files?

Any help with this would be greatly appreciated! Please let me know if I've failed to include any useful information.

Answer: [VCFtools][1] will do this:

> vcftools --recode --vcf your_input.vcf --max-missing (proportion missing above which to filter: [0-1]) --out your_output.vcf

Edit: if you want to use this command per population, you're right: split first.
[1]: https://vcftools.github.io/man_latest.html
Question: Hi!

Fgsea has recently moved to using fgseaMultilevel by default and from what I understand does not use sampling for p-value calculations, but instead an "adaptive multilevel splitting Monte Carlo approach. " (from the documentation). I have been trying to understand some of the new arguments included when running fgsea, particularly **sampleSize**

(`sampleSize` - *The size of a random set of genes which in turn has size = pathwaySize*)

fgseaMultilevel(
pathways,
stats,
sampleSize = 101,
minSize = 1,
maxSize = Inf,
eps = 1e-10,
scoreType = c("std", "pos", "neg"),
nproc = 0,
gseaParam = 1,
BPPARAM = NULL,
nPermSimple = 1000,
absEps = NULL
)

I have been wrestling with the pre-print (https://www.biorxiv.org/content/10.1101/060012v2.full) and all I can gather (might not be accurate) is that this

- has to be an odd number
- it can't be smaller than 3 (that is from their code
https://github.com/ctlab/fgsea/blob/master/R/fgseaMultilevel.R)
- defines the number of gene sets used to calculate probabilities?

But the concept of this `sampleSize` number still remains quite abstract to me (why is it 101 by default), yet it looks very important and is easily changeable (not very foolproof if someone like me plays around with the parameters). And the definition states ...*size = pathwaySize* (confused)?

Basically I am hoping for any input on what does it do in simple terms and when should I change it? The results do vary quite a bit if I try say 11 vs 101 vs 1001 (in terms of the number of significantly enriched pathways). If it's something similar to sampling random gene sets for P-value calculation (which I think is the classical approach) then I am inclined to think this IS something that should potentially be tailored to suit specific needs...

Please help, thanks.

Yours truly,

Desperate For Help

Answer: This is a parameter that control estimation accuracy, somewhat similar to `nperm` in `fgseaSimple`. The higher is the value, the more accurate are the results and the slower (proportionally) it works.

To understand how it's related to accuracy you can use `multilevelError(pval, sampleSize)` function, which tells what would be the estimation error for the given true p-value and used `sampleSize`:

> multilevelError(1e-15, 101)
[1] 1.017545

The result is on log2 scale, so for the above example if the true p-value is 1e-15, `fgseaMultilevel` with the default `sampleSize=101` will give you an estimation from 4\*1e-15 to 0.25\*1e-15 with a 95% probability.

When should you want to change it? Likely never, unlike with `nperm` and `fgseaSimple`, no `sampleSize` value limits how small P-values can be estimated (this can be limited by `eps` argument). The default value of 101 proved to be a good compromise between speed and accuracy. Again, for the given example, I don't know a practical case, when you want to accurately distinguish 1e-15 and 2e-15, while 1e-15 will be distinguishable from 1e-16 and 1e-14. However, if you scale `sampleSize` down to 25 you'll get a 4-fold speed improvement, but the log2 error for 1e-15 will be 2, which means that your CI will be from 16\*1e-15 to 1/16\*1e-15.
Question: I have single-end smallRNAseq data and I have to trim the adapter (adapter sequence: `AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC`). After trimming the adapter using the bbduk command mentioned below, I checked the reads in my data and they still have the overhangs of the adapter as shown below.

The adapter overhangs can be seen below:

cat bbduk_trimmed_small_RNA_001.fastq | grep TCGCAGGGAAATCATCTGATTA

TCGCAGGGAAATCATCTGATTAGATCGGAAGA
TCGCAGGGAAATCATCTGATTAGATCGGAA
TCGCAGGGAAATCATCTGATTAAGATCGGAAGAA
TTCGCAGGGAAATCATCTGATTAGAA

or can be seen here:
https://postimg.cc/image/knuvtyv0d/

This trimming was done using bbduk.sh (bbmap tools: https://jgi.doe.gov/data-and-tools/bbtools/bb-tools-user-guide/bbduk-guide/). Here NEB-SE.fa file has the adapter sequence mentioned above.

bbduk.sh -Xmx1g in=small_RNA_001.fastq out=bbduk_trimmed_small_RNA_001.fastq ref=NEB-SE.fa ktrim=r k=13 mink=6 minlength=18 hdist=0

I could reduce the kmer size `k=4`, but that would risk into trimming the false positives. How can I completely get rid of these adapter sequences from my data without trimming false positives?

Answer: You can't have it both ways. You can either be strict about trimming and risk losing a few reads or have a few reads left that may have a bit of adapter. Your aligner should deal with any extraneous bases while it does the alignment. Have you tried to align this data to see what fractions aligns? With smallRNA you are looking for a specific length (~21-25 bp). If you reads are failing that length criteria after trimming you may want to discard them since they may not be what you are looking for.

BTW: I am curious as to what the original reads in the example you included look like (or is that the full original sequence)?
Question: I'm using the R package 'biomaRt' to retrieve ensembl ID's and descriptions for a list of gene symbols (eg "Sf1", "Rhox7a" etc. . . ). My query consists of 41203 symbols; biomart returns a result set with 30774 records corresponding to the gene symbols recognized by ensembl. The 30774 records returned included **four genes that were not part of the original query**.

My first thought was that the four 'mystery' genes were synonyms for something in my original query. I've since verified that none of the synonyms of these genes are in my query.

I am querying the mouse data set, and using the attribute 'external_gene_name' as my filter column
Code used to query biomaRt

# 'gq': list of unique 'GeneID' submitted as biomart query
gq<-unique(dg$gene)

# attributes used for query
attr<-c("ensembl_gene_id", "external_gene_name", "description",
"ensembl_gene_id_version", "chromosome_name",
"gene_biotype"
)

# Query Submission
mart<-useMart(biomart="ensembl", dataset="mmusculus_gene_ensembl")
result<-getBM(mart=mart,
attributes=attr,
filters='external_gene_name',
values=gq
)

The mystery genes are:

setdiff(result$external_gene_name, gq)
[1] "Trdd2" "Trdv4" "Trdd1" "SPATA24"

Where "gq" is the list of genes submitted to ensembl. None of the above genes, nor any of their synonyms (synonyms recognized by ensembl at least) are in my original query. If anyone is willing to help me troubleshoot, I would be happy to send them the gene list I'm querying with.

Answer: "SPATA24" doesn't look like a normal MGI symbol since it's all in caps, so I wouldn't be suprised if your query contains "Spata24" and the all caps version is retrieved too. Is it possible your gene list include capitalised versions of the 'Trdd1' etc? I don't think BioMart is case senstive and will still retrieve results for them e.g.

getBM(mart=mart,
attributes=attr,
filters='external_gene_name',
values="TRDV4"
)
ensembl_gene_id external_gene_name
1 ENSMUSG00000076867 Trdv4

If that's not it, my advice would be to break your query down into smaller chunks and submit this independently, to try and narrow down where the unexpected entries are being introduced. Happy to try and identify if it's a problem in **biomaRt**, email address is on the
**biomaRt** landing page (https://bioconductor.org/packages/biomaRt/)
Question: I'm trying to automatically compare BAM files being output by bowtie2 (for a continuous integration system). Currently using `bamUtil diff` which works well, except I'm noticing sometimes 2 lines in the BAM file will have the same name/start position, but represent different reads. For whatever reason, bowtie2 sometimes reverses the order of the lines, so bamUtil's algorithm calls a mismatch. I see several potential solutions to this problem:

1. Get samtools to sort by another dimension in addition to position/name. Pretty much anything should do. Can this be done? I don't see any options in the man page.

2. Use another comparision tool than bamUtil which doesn't have this problem. Honestly a stand-alone program for diffing would be best for me anyway and I'm considering writing one.

3. Find some way to get bowtie2 to be more deterministic in its ordering. Any ideas here?

Answer: I kinda feel like an idiot right now. Upon further reflection on my task, I realized that in this particular case there's no reason to limit myself to samtools' method of sorting. I can simply convert my two BAM files to SAM (without the headers), then feed them to GNU coreutils sort, then diff them normally. Seems to be working. Here's a bash script that does it:

#!/bin/bash
# bamdiff.sh

f1=${1}
f2=${2}

f1TempHeaderlessSAM=.bamdiff_f1TempHeaderlessSAM
f2TempHeaderlessSAM=.bamdiff_f2TempHeaderlessSAM

samtools view ${f1} | sort > ${f1TempHeaderlessSAM}
samtools view ${f2} | sort > ${f2TempHeaderlessSAM}

diff ${f1TempHeaderlessSAM} ${f2TempHeaderlessSAM}

diffExitCode=$?

rm ${f1TempHeaderlessSAM}
rm ${f2TempHeaderlessSAM}

exit ${diffExitCode}
Question: Dear everyone,

I am doing analysis from GEO database, I download the soft files, extract the expression matrix of probes-to-sample but I don't know where to find clinical data for these sample for further survival analysis. Could you please guide me how to do this.

For example this GSE data: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=gse10780

Thank you!

Answer: Reproducibility is huge problem and most distribution systems like GEO are not well suited for distributing data other than a few simple and known types like fastq, fasta etc.

Here is what you can do:

1. Check the supplementary data in the paper - these typically are notoriously obtuse and convoluted to deal with. Often are not described in sufficient detail hence you will need to download and open each to see what is inside them.
2. Your actual best bet is to contact the authors and ask for this data (or have them point out where this data is).
Question: Hi,

I have a BAM file which has chr id's as `NC_00000*.` I did sorting using samtools sort function.

I wanted to remove duplicates, so I'm using MarkDuplicates.jar from Picard tools to get the job done. But it gives me the following error:

Exception in thread "main" net.sf.picard.PicardException: 13_0501.sorted.bam is not coordinate sorted.
at net.sf.picard.sam.MarkDuplicates.buildSortedReadEndLists(MarkDuplicates.java:273)
at net.sf.picard.sam.MarkDuplicates.doWork(MarkDuplicates.java:117)
at net.sf.picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:158)
at net.sf.picard.sam.MarkDuplicates.main(MarkDuplicates.java:101)

But I think my bam file is sorted. This is the header of my bam file.

@HD VN:1.0 SO:unsorted
@SQ SN:NC_000001 LN:249250621
@SQ SN:NC_000002 LN:243199373
@SQ SN:NC_000003 LN:198022430
@SQ SN:NC_000004 LN:191154276
@SQ SN:NC_000005 LN:180915260
@SQ SN:NC_000006 LN:171115067
@SQ SN:NC_000007 LN:159138663
@SQ SN:NC_000008 LN:146364022
@SQ SN:NC_000009 LN:141213431
@SQ SN:NC_000010 LN:135534747
@SQ SN:NC_000011 LN:135006516
@SQ SN:NC_000012 LN:133851895
@SQ SN:NC_000013 LN:115169878
@SQ SN:NC_000014 LN:107349540
@SQ SN:NC_000015 LN:102531392
@SQ SN:NC_000016 LN:90354753
@SQ SN:NC_000017 LN:81195210
@SQ SN:NC_000018 LN:78077248
@SQ SN:NC_000019 LN:59128983
@SQ SN:NC_000020 LN:63025520
@SQ SN:NC_000021 LN:48129895
@SQ SN:NC_000022 LN:51304566
@SQ SN:NC_000023 LN:155270560
@SQ SN:NC_000024 LN:59373566
@PG ID:0 PN:clcgenomicswb VN:7.0

The header seems to suggest that it's unsorted. That's what is bothering me.

Answer: your SAM Header starts with:

@HD VN:1.0 SO:unsorted

And Picard check that this 'SO:' is set to 'coordinate'. See https://github.com/nh13/picard/blob/master/src/java/net/sf/picard/sam/MarkDuplicates.java

(....)
if (!ASSUME_SORTED && header.getSortOrder() != SortOrder.coordinate) {
throw new PicardException("Input file " + f.getAbsolutePath() + " is not coordinate sorted.");
}

to fix this, use the option

ASSUME_SORTED=true
Question: I came across a lot of literature that uses 30x reads as a good coverage depth for NGS. But then I read up on Broad's NGS intro where they used "high-pass: 30x; low-pass: 4x; exon capture:150x". Not sure when the intro note was written, so no idea how relevant those numbers are to HiSeq 2000 machine.

Q1. how do they come up with those numbers? (anything to do with their Bayesian model-based SNP calling ?)

Q2. how are those numbers related to detecting rare SNPs in a small pool of sample? Say, I want to find out a SNP whose population-based occurrence is between 0.1-0.01%. I only have a handful of individual samples (say 20). Should I adjust the coverage depth (say, from 30x to 100x) in order to detect with high confidence the SNP of interest?

Some digging around in google and pubmed does not really help me find the answers. Any help/reference will be greatly appreciated. Thank you.

Answer: This is a hard question to answer definitively. It really depends on what you are trying to do and especially what kind of samples you are profiling. What is the expected purity of your sample? Is there contamination or heterogeneity expected? If the answer to those last two questions is yes, then you will want more coverage to detect variants at lower frequency due to impurity or because they come from minor (but possibly still interesting) subpopulations of your sample. We tend to target 200-250x for tumor studies for the reasons described above. But, I don't think there is a really solid statistical basis for that choice. The cost-benefit analysis is probably too complex for a definitive answer so practical considerations such as cost tend to take greater precedence.

If you are looking for a very rare SNP in a population of individuals (rather than a population of cells in a single sample) then I don't think extra coverage is nearly as important as sampling more individuals. If a SNP is only found in 0.01% of the population then it is highly unlikely that you will detect it in a sample of 20 individuals no matter how deeply you sequence. Chances are the SNP just isn't there. Am I missing something in your question?

Here are other related BioStar questions:

What depth of coverage for whole genome SNP calling?

What coverage is required for ChIP-Seq experiments?

How much coverage do we need for an RNA-seq experiment?
Question: hi,

I have a bed file like below

```
FLP1 252 1523
RAF1 3271 3816
REP1 1887 3008
REP2 5308 6198
```

I was reading the file in R

```
bed <- as.data.frame(read.table("file.bed",header = FALSE, sep="\t",stringsAsFactors=FALSE))
Warning messages:
1: In scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings, :
EOF within quoted string
2: In scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings, :
number of items read is not a multiple of the number of columns
```

but when I read file as csv no error anymore

how I can read this file as table please?

thank you

Answer: Maybe you have a weird hidden character in there somewhere. try adding

quote = ""

to your initial `read.table` call, like:

bed <- as.data.frame(read.table("file.bed",header = FALSE, sep="\t",stringsAsFactors=FALSE, quote=""))