cannin/biostars_qa huggingface dataset:
Question: Dear all,

I am trying to install [Primer3][1] on an Ubuntu machine. Following the instructions, I moved into primer3-/src and typed `make all` but i got:

/usr/local/lib/primer3/src$ make all
g++ -c -g -Wall -D__USE_FIXED_PROTOTYPES__ -O2 primer3_boulder_main.c
g++ -c -g -Wall -D__USE_FIXED_PROTOTYPES__ -O2 -o format_output.o format_output.c
g++ -c -g -Wall -D__USE_FIXED_PROTOTYPES__ -O2 -o read_boulder.o read_boulder.c
g++ -c -g -Wall -D__USE_FIXED_PROTOTYPES__ -O2 -o print_boulder.o print_boulder.c
g++ -c -g -Wall -D__USE_FIXED_PROTOTYPES__ -O2 -Wno-deprecated -o libprimer3.o libprimer3.c
g++ -c -g -Wall -D__USE_FIXED_PROTOTYPES__ -O2 -o p3_seq_lib.o p3_seq_lib.c
ar rv libprimer3.a libprimer3.o p3_seq_lib.o
ar: creating libprimer3.a
a - libprimer3.o
a - p3_seq_lib.o
ranlib libprimer3.a
g++ -c -g -Wall -D__USE_FIXED_PROTOTYPES__ -O2 -o dpal_primer.o dpal.c
dpal.c: In function ‘void print_align(const unsigned char*, const unsigned char*, int (*)[1600][3], int, int, const dpal_args*)’:
dpal.c:1036:5: warning: this ‘for’ clause does not guard... [-Wmisleading-indentation]
for(i=j;i

Answer: I finally was able to install it following the instructions in the github page here https://github.com/primer3-org/primer3

Installing
sudo apt-get install -y build-essential g++ cmake git-all

git clone https://github.com/primer3-org/primer3.git primer3

cd primer3/src

make

make test
Question: I have a bed file called `my.bed` with CHROM, START, and END Position. Can someone please explain me how I can use bcftools or command to extract the regions from `myvcf.vcf` file?

Answer: my preferred option when dealing with vcf files is [bcftools][1], which requires vcf indexing:

tabix -p vcf my.vcf
bcftools view -R my.bed my.vcf.gz

another perfectly valid alternative, as Devon has just pointed out, would be [bedtools intersect][2]:

bedtools intersect -a my.vcf -b my.bed

[1]: http://samtools.github.io/bcftools/bcftools.html#view
[2]: http://bedtools.readthedocs.org/en/latest/content/tools/intersect.html
Question: Hello,

I am trying to get haplotypes for one locus on chromosome19 and for that, I need phased data. However, prior to the phasing, I am supposed to split my PLINK files into separated files by chromosome. I found this script online but it doesn't work (see below). Can you please give me advice on how to do this? I am rather a newbie in the data processing so I am sorry if this is a too basic question.

Thank you!

#!/usr/bin/perl

# This script takes as input the base filename of binary pedfiles (*.bed,
# *.bim, *.fam) and a base output filename and splits up a dataset by
# chromosome. Useful for imputing to 1000 genomes.

chomp(my $pwd = `pwd`); my $help = "\nUsage: $0 \n\n"; die $help if @ARGV!=2;

$infile_base=$ARGV[0]; #base filename of inputs $outfile_base=$ARGV[1]; #base filename of outputs $plink_exec="plink
--nonfounders --allow-no-sex --noweb"; $chr=22; #last chromosome to write out

for (1..$chr) {	print "Processing chromosome $_\n";	`$plink_exec
--bfile $infile_base --chr $_ --make-bed --out ${outfile_base}$_;` }

Answer: I think the correct format is (in case some codes are commented out):

```
#!/usr/bin/perl

# This script takes as input the base filename of binary pedfiles (*.bed,
# *.bim, *.fam) and a base output filename and splits up a dataset by
# chromosome. Useful for imputing to 1000 genomes.

chomp( my $pwd = `pwd` );
my $help = "\nUsage: $0 \n\n";
die $help if @ARGV != 2;

$infile_base = $ARGV[0]; #base filename of inputs
$outfile_base = $ARGV[1]; #base filename of outputs
$plink_exec = "plink --nonfounders --allow-no-sex --noweb";
$chr = 22; #last chromosome to write out

for ( 1 .. $chr ) {
print "Processing chromosome $_\n";
`$plink_exec --bfile $infile_base --chr $_ --make-bed --out ${outfile_base}$_;`;
}
```
Question: Hi everyone,

My lab is using gene expression data generated by Illumina Human HT-12 v3 Expression Beadchips. As advertised by the company, this products has 48000+ probes for 25000 genes. I have never used expression data before and would like to cluster genes based on their expression. The data has already been normalized and corrected for batch effects.

The current file format is:

ProbeID Sample1 Sample2

I would like to get the following format:

GeneID Sample1 Sample2

It seems that some genes have more probes than others. Moreover, there can be multiple transcripts for a given gene. I was wondering if someone could please give me a general idea about getting the desired format.

Thank you for your time.

Answer: Hi,

Its easier to do this in R. All you need is to convert ProbeID into the Gene name to which it is mapped.

```
> probeID=c("ILMN_1690170", "ILMN_2410826", "ILMN_1675640", "ILMN_1801246",
"ILMN_1658247", "ILMN_1740938", "ILMN_1657871", "ILMN_1769520",
"ILMN_1778401")
> library("illuminaHumanv4.db") #Get this library if you don't have

> data.frame(Gene=unlist(mget(x = probeID,envir = illuminaHumanv4SYMBOL)))
Gene
ILMN_1690170 CRABP2
ILMN_2410826 OAS1
ILMN_1675640 OAS1
ILMN_1801246 IFITM1
ILMN_1658247 OAS1
ILMN_1740938 APOE
ILMN_1657871 RSAD2
ILMN_1769520 UBE2L6
ILMN_1778401 HLA-B
```
Question: Dear All,

I have a big embl file that has over 1000 IDs, and I want to split it based on IDs. Each ID should be a separate file containing all information of related ID, and should have the ID name.

How to do that?

Answer: @OP: please post first embl record for better parsing.

Note: Before proceeding, create a test directory, copy the original embl (test.embl below) and run the script and check the output for random files.

OP embl ID line:

> ID scaffold00001; SV 1; linear; unassigned DNA; STD; UNC; 6279 BP.

based on OP ID line,

$ grep -i "ID" test.embl | awk '{gsub(";",""); print $2}' | parallel "awk '/{}/,/\/\//' test.embl > {}.embl"

Above script should create multiple embl files with ID as file name.

Example embl file from post https://www.biostars.org/p/147163/ (close to OP embl IDs):

$ cat test.embl
ID comp0_c0_seq1; SV 1; linear; unassigned DNA; STD; UNC; 205 BP.
XX
DE len=205 path=[1:0-135 1445:136-204]
XX
SQ Sequence 205 BP; 64 A; 54 C; 31 G; 56 T; 0 other;
GTATTGAACT GCAGAGCATT AAATGCTGCA ACTCAGTGCT TAGAATTCAT TAGATTCAGA 60
GCAACGAACC CTAAATACTG AGCTGTCCCA TTAAATACTC TGCAGTTCAA TACTTAGCAT 120
TCACCATTAA ACATAACACT TCCCGAGTTT CCACCATCCA TAAACAGCAG GCATTGTAAC 180
CTGTAGGCTC TCTCCACGGT TACCT 205
//
ID comp0_c0_seq2; SV 1; linear; unassigned DNA; STD; UNC; 205 BP.
XX
DE len=205 path=[4094:0-135 1445:136-204]
XX
SQ Sequence 205 BP; 59 A; 50 C; 35 G; 61 T; 0 other;
AGAGTATTAA ATGTTGCAGT TCAGTGCTTA AAATTTATTG GATTCAGAGA ATCTTCAAAT 60
TCAACGGACC CTAAACACTG AGCTGTCGCA TTAAATGCTC TGCAGTTCAA TGCTTAGCTT 120
TCACCATTAA GCATAGCACT TCCCGAGTTT CCACCATCCA TAAACAGCAG GCATTGTAAC 180
CTGTAGGCTC TCTCCACGGT TACCT 205
//
ID comp1_c0_seq1; SV 1; linear; unassigned DNA; STD; UNC; 244 BP.
XX
DE len=244 path=[3:0-88 875:89-243]
XX
SQ Sequence 244 BP; 71 A; 51 C; 63 G; 59 T; 0 other;
GCAGAATTTA AGGCTATGAA TCAGGAGGTT CATAATTCCT TAAGGAGGGG AGTATGATGC 60
GGAGCATCCA CGCTCACCTC CACTCCACCG CATTGTCTTC GAGCTGTGAC AGCCAGCGCA 120
TAATATTCAA GAGCTATTGA CAGGTGTTGA AACGCGGCAG CCTTGCATAC TATTGAAGGA 180
CCACGTTTCA TTATTGTGAT CTATAAGAAG ACAGCTGATG CGATCATGAG GAAGGAAGAA 240
GGCT 244
//

output:

$ ls *.embl
comp0_c0_seq1.embl comp0_c0_seq2.embl	comp1_c0_seq1.embl test.embl
.

$ more comp0_c0_seq1.embl
ID comp0_c0_seq1; SV 1; linear; unassigned DNA; STD; UNC; 205 BP.
XX
DE len=205 path=[1:0-135 1445:136-204]
XX
SQ Sequence 205 BP; 64 A; 54 C; 31 G; 56 T; 0 other;
GTATTGAACT GCAGAGCATT AAATGCTGCA ACTCAGTGCT TAGAATTCAT TAGATTCAGA 60
GCAACGAACC CTAAATACTG AGCTGTCCCA TTAAATACTC TGCAGTTCAA TACTTAGCAT 120
TCACCATTAA ACATAACACT TCCCGAGTTT CCACCATCCA TAAACAGCAG GCATTGTAAC 180
CTGTAGGCTC TCTCCACGGT TACCT 205
//
$ more comp1_c0_seq1.embl
ID comp1_c0_seq1; SV 1; linear; unassigned DNA; STD; UNC; 244 BP.
XX
DE len=244 path=[3:0-88 875:89-243]
XX
SQ Sequence 244 BP; 71 A; 51 C; 63 G; 59 T; 0 other;
GCAGAATTTA AGGCTATGAA TCAGGAGGTT CATAATTCCT TAAGGAGGGG AGTATGATGC 60
GGAGCATCCA CGCTCACCTC CACTCCACCG CATTGTCTTC GAGCTGTGAC AGCCAGCGCA 120
TAATATTCAA GAGCTATTGA CAGGTGTTGA AACGCGGCAG CCTTGCATAC TATTGAAGGA 180
CCACGTTTCA TTATTGTGAT CTATAAGAAG ACAGCTGATG CGATCATGAG GAAGGAAGAA 240
GGCT 244
//
Question: Hello, I am wondering, what tool/R-package can be used to draw this type of tree (preferably, programmatically)?

![enter image description here][1]

The image is taken from: https://www.ncbi.nlm.nih.gov/pubmed/26773003

[1]: http://ivanya.com/tmp/2016/1024.tree.png

Answer: Now it is supported by [ggtree](https://guangchuangyu.github.io/2016/11/align-genomic-features-with-phylogenetic-tree/).
Question: In Biostars forum for differential gene expression analysis by LIMMA I found that most people suggest to set a threshold value of Log2FC > or = 2 to filter DEGs. When I analyzed the GEO dataset "GSE90594" - Agilent-039494 SurePrint G3 Human GE v2 8x60K Microarray 039381 platform for DEGs by LIMMA package, I found that very few genes only have a log 2FC value above 1.

But in the original paper "Study of gene expression alteration in male androgenetic alopecia: evidence of predominant molecular signalling pathways" (http://doi.wiley.com/10.1111/bjd.15577), the authors reported that they obtained 325 UP and 390 DOWN regulated DEGs. They have followed a different approach for obtaining DEGs and reported the DEGs with Fold change (FC) values. Most of the UP regulated genes are reported with the FC value 1.5 to 2 range (log2FC value 0.58 to 1.3 approx). Similarly for down regulated genes FC values of 0.5 to 0.8 (log2Fc value -1 to -0.3) are considered.

Is this right ? can we consider log2FC value of above 0.6 for UP and below -0.6 for DOWN DEGs. Please someone clarify this irrespective of p and FDR values. I consider a q value of 0.05 for DEG selection.

Answer: The fold change will tell you something of the size of the effect of differential expression, so it's more about biology than about statistics. If you are looking for genes with big changes, you'll pick a higher cut-off. But if a gene is significantly differentially expressed then it's already worth looking into it: subtle differences in gene expression can already have a substantial impact. It probably makes a big difference in which gene is affected, too. Some genes are more dosage-sensitive than others.
Question: I have RNA-seq data with 2 conditions and 3 replicates per conditions.

I ran the [New Tuxedo pipeline][1] and also created some [read count tables][2] with prepDE.

I analysed differentially expressed genes with `Ballgown` and `DESeq2`.

With a treshold of 1 log2FoldChange and 0.01 padj in `DESeq2`: 14400 /32000 (45%) of DE genes

With a treshold of 0.01 pval in `Ballgown`: 3678/32000 (of DE genes), even with no fold change treshold, the number of DE genes is (very) lower. In ballgown, what is the difference between qval and pval ? Which one corresponds to padj in DESeq2 ?

I expect many DE genes as conditions are very different biogically (testis vs ovary, same species).

Why do I have a so large difference between softwares?

[1]: https://www.nature.com/articles/nprot.2016.095
[2]: http://ccb.jhu.edu/software/stringtie/index.shtml?t=manual

Answer: *pval* is the nominal p-value. *qval* is the adjusted p-value, which are also known as q-values (not many people know this).

Ballgown may be using FPKM data when conducting the differential expression analysis. FPKM is not suitable for this purpose. Please confirm the type of normalisation that you used in Ballgown.

When you ran DESeq2, did you use the `lfcShrink()` function? - see the piece of code that I posted here: https://www.biostars.org/p/324128/#324305
Question: Hi all, I am new to the BioInformatics, and quite a beginner in programming languages. Can anyone suggest me some sources where I can at least learn 50% of the scRNA seq data analysis? I am familiar with C language, and I know a little bit of molecular biology too.

Answer: Single-cell data are rather unpleasant as a beginner's topic due to the noisy and sparse nature of these data. Maybe better first analyze some bulk RNA-seq data to get familiar with R (see [here][1]), and then dive into the documentation of `Seurat` which is the jack-of-all-trades in terms of scRNA-seq analysis. For lowlevel processing `alevin` is a good choice.

[1]: https://www.bioconductor.org/packages/devel/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html
Question: Hello,

I don't know why I am getting some errors during my analysis

I uploaded an example of my data in

https://gist.github.com/anonymous/2c69ab500bfa94d0268a

I use the following command in R to load my data

data <- read.delim("path to your file /example.txt", header=FALSE)

however, in summary or head or other commands I look at the data, it seems alright but I cannot analysis since it gives error like `all numeric variables`. For example if you want to get the range of the example data, you will get such error.

How normally do you import, load a microarray data (with **txt format**) (each row represents a prob and each column a sample)?

Thanks

Answer: Works perfectly for me:

```
> dat <- read.delim(file="Downloads/gist2c69ab500bfa94d0268a-ac4cd3d5b0d0764c2faae0e3fb0db8a39d75bb22/example.txt", row.names=1)
# your mistake was to set header=FALSE, and to omit
# row.names=1

> head(dat)
M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 M12
200645_at 0.0446 0.0744 -0.0340 0.0173 0.2280 0.0070 -0.0250 0.0644 -0.0253 -0.1230 -0.6251 0.0210
200690_at -0.0165 0.1121 -0.0959 0.0000 -0.4595 -0.0282 -0.1617 -0.0482 -0.2611 0.0223 -0.6129 0.1961
200691_s_at 0.0554 -0.0689 -0.0852 0.0702 0.0823 0.0361 -0.0306 -0.0076 -0.0340 -0.0198 -0.1823 -0.0681
200692_s_at 0.0000 -0.0505 -0.0508 -0.0159 -0.3041 -0.0684 -0.0644 -0.0175 0.0503 0.0546 -0.2141 -0.0216
200693_at 0.0608 0.0601 0.0115 0.0744 -0.0232 -0.1095 -0.0416 -0.0499 -0.0515 0.0303 -0.1153 0.0824
200694_s_at 0.0424 0.0957 0.0758 -0.0387 -0.0517 -0.0207 0.0328 -0.1392 0.0140 -0.1476 0.1382 0.0113
M13 M14 M15 M16 M17 M18
200645_at 0.1095 0.1527 0.0261 -0.2107 -0.0196 -0.2316
200690_at 0.2119 0.0122 -0.5495 0.1518 -0.2409 0.1610
200691_s_at 0.1219 -0.1615 -0.0729 -0.0696 0.0042 0.1239
200692_s_at 0.0440 -0.0811 0.0964 0.0211 -0.0325 0.1810
200693_at -0.0036 0.0575 0.0427 0.1104 -0.0216 0.0278
200694_s_at 0.2247 0.1489 0.0196 0.0883 -0.1848 0.1989

> range(dat)
[1] -20.091 25.652
```
Question: What exactly are RMA-units I get after normalizing affymetrix data?

These are relative expression log2 units. But relative to what? What is RMA = 0?

Answer: See this link to a presentation:

http://www.ub.edu/stat/docencia/bioinformatica/microarrays/ADM/slides/2_PreprocessingMicroarrayData-2-Preprocessing%20and%20Normalization.pdf

And also this link to some discussion on the topic:

https://www.researchgate.net/post/What_exactly_is_the_unit_of_measure_in_microarray_experiments

" Michael B Black · ScitoVation
The raw data is effectively unitless as it is simply relative signal intensity. Flourescence is measured by a photo multiplier tube or charge-coupled device and signal scaled across the range of detection for the platform. The starting data is thus simple relative signal values without units, based effectively on simple counted photons by the PMT or CCD.
"

" Prafull Kumar Singh · University Medical Center Freiburg
Micro array chip contain multiple probes for each gene. After completion of the experiment a .dat is generated which is a image showing the intensity of each probe in a digital format. Now it is converted into .cell file which contains the real probe intensity numerical values.The higher the probe intensity higher is the expression of the gene. But its meaningless until you compare it with the intensity of same probe in control sample. These probe intensity values are preprocessed (Background correction, normalization, Perfect match correction and Summarization) using diff algorithms (RMA, Mas5). The normalized probe intensity values are further used for calculation of relative gene expression i.e treatment vs control to get the fold change. The positive and negative values indicates up-regulated and down-regulated genes in treatment group when compared to control."

etc...

An initial experiment reference:

http://biostatistics.oxfordjournals.org/content/4/2/249.full.pdf
Question: Hi,

This appears to be a simple problem that I am unable to solve. I have some data that looks like this:

CHROM POS REF ALT TYPE AF
chr1 1 A T MISSENSE 0.23
chr2 1 A T,G MISSENSE 0.17, 0.09

The above is dummy meaningless data, but it is representative of the problem at hand.

I'd like to `separate_rows` such that the `ALT` and `AF` are separated in a couples manner. Running `separate_rows` on the 2 columns would give me 4 rows, not 2. I'd like my output to be:

CHROM POS REF ALT TYPE AF
chr1 1 A T MISSENSE 0.23
chr2 1 A T MISSENSE 0.17
chr2 1 A G MISSENSE 0.09

Is there any way I can conserve this combination while separating the values out? I am really far out from the VCF to go back and split multi-allelics.

Answer: You could use a Python script to do this easily:

#!/usr/bin/env python

import sys

headers = None
idx = 0
for line in sys.stdin:
elems = line.rstrip().split('\t')
if idx == 0:
headers = elems
sys.stdout.write(line)
else:
items = {x:y for x,y in zip(headers, elems)}
alleles = items['ALT'].split(',')
afs = items['AF'].split(',')
for ai in range(len(alleles)):
items['ALT'] = alleles[ai]
items['AF'] = afs[ai]
sys.stdout.write('{}\n'.format('\t'.join([items[x] for x in headers])))
idx += 1

For example:

$ ./split.py < variants.txt
CHROM	POS	REF	ALT	TYPE	AF
chr1	1	A	T	MISSENSE	0.23
chr2	1	A	T	MISSENSE	0.17
chr2	1	A	G	MISSENSE	0.09

Write it out to a file and bring that back into R:

$ ./split.py < variants.txt > variants.split.txt
Question: Hello everyone,

I am really new in this forum, so I hope you can help me in resolving my problem.
I am actually working on a project, where I have to annotate some sequences.these sequences are a result of rna_seq , and I have them in bam format. Do you know some easy bioinformatics tools that can help me to annnotate that kind of sequences in order to get GFF files ?? So the input should be in bam files and the output in GFF.
Just for information, i need a tool that can be executed in linux, and that can be looped for so many sequences for different species.

Thank you for your help

Answer: OK, so you want to annotate (predict genes) on a genomic sequence making use of RNAseq info.

The question is easily asked but I'm afraid you might be underestimating how a complex matter this can be. But is absolute a good and correct approach. Problem is that this often requires quite a bit of work and knowledge. There are (or you could) run most tools as a black box but then your results will be of lesser quality.

You might consider tools as: Maker, Braker, Augustus, ... (most of these use RNAseq as an evidence information but also require other inputs ). There will be others that might work with RNAseq only but none come to mind right now.
Question: Hi,

I have been trying to use multiple command i can find online to spilt the columns into different files but couldnt.

I have an 10 columns excel file ( first 2 col is paramaters, the remaining 8 col is the data). I would like to split it into 8 separate files which all 8 contains first 2 col paramaters and 1 data col.

May i know does anyone have any idea to get around this ?

Thank you.

Answer: Here is the code in R: Make sure that all for the files (excel files with .xls or .xlsx) are in the same format as in link furnished above (2nd column is mz, 5 is RT and 11th column onward samples). Code would create one excel file per sample and name of the excel file would be sample_"sample name".xls. It will have three columns: mz, RT and sample name.

library(readxl)
library(WriteXLS)
df=data.frame(read_xls("example.xls"))

for (i in 11:ncol(df)){
temp_df=df[,c(2,5,i)]
WriteXLS(temp_df, ExcelFileName = paste0("sample_",names(df)[i],".xls"))
}

if you have multiple xls files in the folder (and only xls files of interest):

setwd("~/Desktop/test") # change this to directory of interest
fn=list.files(pattern = "\\.xls") # lists files with .xls extension
library(readxl)
library(WriteXLS)

for (i in fn){
df=data.frame(read_xls(i))
for (j in 11:ncol(df)){
temp_df=df[,c(2,5,j)]
WriteXLS(temp_df, ExcelFileName = paste0(sub('\\.xls$', '', i) ,"_","sample_",names(df)[j],".xls"))
}
}

output files will have "xlsfilename_samplename.xls"
Question: Hello, I am looking at this heatmap and I do not understand why some of the tumours are grouped with the controls, seems as if the heatmap is 'moved to the left':
Top of the heatmap:
![Top of the heatmap][1]

(it is a long heatmap so i upload just a section)

Bottom of the heatmap:
![Bottom of the heatmap][2]

I tried to check on the phenotype of the samples but I found no correlations of the phenotype with this grouping.

There seem to be two tumour subgroups: one that is the biggest, and the small group of tumors that seem to be grouped with the controls. I did a t-test on the mean beta values for all the cpgs between these two groups, and it turns out there is significance among the means of these two tumour groups.
I am afraid when we publish a heatmap similar to this one, we will have trouble explaining this phenomenon. Any ideas or any opinions on this? Thank you!

[1]: https://i.imgur.com/v3oAPra.png
[2]: https://i.imgur.com/NgTGrB3.png

Answer: There can be many reasons, some not really relating to bioinformatics:

1. these tumours genuinely exhibit the 'normal-like' methylation
profile over these probes
2. these tumours have normal cell contamination

Some informatics reasons:

- your coding is incorrect and you have incorrectly assigned a normal
sample as a tumour
- you have scaled your data incorrectly
- you should consider the distance and linkage metric that you're using

By the way, for methylation, you could probably also perform the Wilcoxon Signed Rank test on the matched T-N pairs. If just a regular t-test, I would at least use the Mann Whitney (non-parametric) test.

You must also have an extra checkpoint: you should obtain the difference in mean β value between tumour and normal, i.e.:

difference in mean = mean β (tumour) - mean β (normal)

Then, use that as an extra cut-off in addition to the p-value.

Kevin
Question: How can we get total number of SNPs in plink association output file as "qassoc" and "adjusted file" which has less than 0.05
P values?

Answer: From the plink manuals the pvalue column for [qassoc file][1] is the 9th. So we can use `awk` to keep the 1st row header, then filter on 9th column value, see example:

awk '(NR==1) || ($9 < 0.05) ' myfile.qassoc > myfile_subset.qassoc

To get the number of rows:

awk '$9 < 0.05' myfile.qassoc | wc -l

[1]: https://www.cog-genomics.org/plink2/formats#qassoc
Question: I am trying to convert a bigwig file to a wig file to run through sitepro. I tried using the ucsc tool bigWigToWig and I get a file that looks like this:

```
#bedGraph section chr1:0-3079820
chr1 0 3001030 0
chr1 3001030 3001060 2.41
chr1 3001060 3001230 4.83
chr1 3001230 3001270 2.41
chr1 3001270 3002330 0
chr1 3002330 3002380 4.83
chr1 3002380 3002540 9.65
chr1 3002540 3002590 4.83
chr1 3002590 3003210 0
```

However, sitepro is choking on this wig file. I normally use wig files generated by MACS, which have a much different formatting:

```
track type=wiggle_0 name="macs_output_treat_all" description="Extended tag pileup from MACS version 1.4.2 20120305 for every 10 bp"
variableStep chrom=chr1 span=10
3001071 1
3001081 2
3001091 2
3001101 2
3001111 2
3001121 2
3001131 2
3001141 2
```

Why is the ucsc bigWigToWig tool giving me a bedgraph file rather than a wig file?

Answer: This guy experienced the same problem, and has an answer (and a script to go from bedGraph to wig) on his site:

http://sebastienvigneau.wordpress.com/2014/01/10/bigwig-to-bedgraph-to-wig/