cannin/biostars_qa huggingface dataset:
Question: Hi,

I've got a problem with specific bam files, I've downloaded only a slice of them (one chromosome), this results in some `MATE_NOT_FOUND` error when checking my bams with:

gatk ValidateSamFile -I input.bam -M SUMMARY

Giving the error:

```
Error Type Count
ERROR:MATE_NOT_FOUND 1408
ERROR:MISMATCH_FLAG_MATE_NEG_STRAND 18551
ERROR:MISMATCH_FLAG_MATE_UNMAPPED 733
```

I have tried two solutions, both decreases the number of `MATE_NOT_FOUND` but there are still a few remaining, I don't understand how it is possible.

**1\. I've tried with samtools**:

samtools view input.bam -f 0x1 -f 0x2 -b -o input.fixed.bam

It should only keep "read paired (0x1)" and "read mapped in proper pair (0x2)".

Now I have this output:

```
Error Type Count
ERROR:MATE_NOT_FOUND 81
ERROR:MISMATCH_FLAG_MATE_UNMAPPED 524
```

**2\. I've tried with PrintReads (GATK)**:

gatk PrintReads -I input.bam -o input.fixed.bam --read-filter PairedReadFilter

It should only keep reads that are properly paired:

Now I have this output:

```
Error Type Count
ERROR:MATE_NOT_FOUND 81
ERROR:MISMATCH_FLAG_MATE_UNMAPPED 524
```

Why they are not all removed?

Answer: I've tried the solution proposed by @genomax but it didn't work for me.

After many attempts with the bams I finally solved the problem with a more radical solution by reverting the bams to fastq

```
samtools sort -n -O BAM -o sorted.bam input.bam
samtools fastq -1 output_1.fastq.gz -2 output_2.fastq.gz -0 /dev/null -s /dev/null -n sorted.bam
```

Then realign with `bwa mem`.

I think it is not the optimal solution as data loss could occur (BQSR data?) but at least it works.
Question: I am trying to understand the information present in PhastCons elements bed files from USCS Table Browser. Following the information I found [here][1] I managed to get a description of the column names, but not a description of what they are measuring. I checked other BioStars posts, and apparently, the PhastCons have scores that range from 0-100, but this is not true according to the description from Table Browser (see below):

![enter image description here][2]

I also checked the distribution of the the "score" column, and on the file that is available to download, the smallest score is 117. So I am assuming the phastConsElements46wayPrimates bed file only contains the conserved regions. Is this correct? Also, where can I get a description of what the LOD score is measuring here.

[1]: http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/?C=M;O=A
[2]: /media/images/1aae92a6-1ed5-46f1-bb9b-01cb4d6f

Answer: Hello,

We offer more information in the track description (https://genome.ucsc.edu/cgi-bin/hgTrackUi?db=hg19&c=chrX&g=cons46way).

There's 4 components to our conservations tracks:

- A multiz alignment of the sequences
- Basewise conservation as predicted by phyloP - this is displayed as a wiggle track
- Element conservation as predicted by phastCons - this is also displayed as a wiggle, but less resolution than the phyloP
- Conserved elements as predicted by phastCons with the --viterbi option (I believe this is now --most-conserved, it used to be viterbi: http://compgen.cshl.edu/phast/phastCons-tutorial.php) - this is discrete elements in bed format, which is the file you are looking at

Below is the description of that specific track. lod stands for "log-odds score", and the score is ultimately transformed to be between 0-1000 to optimize the Genome Browser's display.

> The conserved elements were predicted by running phastCons with the
> --viterbi option. The predicted elements are segments of the alignment that are likely to have been "generated" by the conserved state of the
> phylo-HMM. Each element is assigned a log-odds score equal to its log
> probability under the conserved model minus its log probability under
> the non-conserved model. The "score" field associated with this track
> contains transformed log-odds scores, taking values between 0 and
> 1000. (The scores are transformed using a monotonic function of the form a * log(x) + b.) The raw log odds scores are retained in the
> "name" field and can be seen on the details page or in the browser
> when the track's display mode is set to "pack" or "full".

If you have any follow up questions, our public help desk can always be reached at genome@soe.ucsc.edu. You may also send questions to genome-www@soe.ucsc.edu if they contain sensitive data. For any Genome Browser questions on Biostars, the UCSC tag is the best way to ensure visibility by the team.
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: Filter first by adjusted p-value and only then by a fold-change that makes sense for your experiment. This might be a fold-change of 1 or 0.5 or 2 or something else, there's no one-size-fits-all solution to that. In some cases you may not filter by fold-change at all.
Question: I have two samples with me. they are CEL files. image file form microarray studies. I would like to analyse and determine the differentially expressed genes . Now my question is anybody please tel me how can I analyse the microarray CEL files?

And is there any free software available for this kind of analysis?

Thank you

Answer: If you have only two samples without replication it is not possible to perform a proper DE analysis using statistical tests. The only thing you can do is rank the genes by absolute differences between those two samples, extracting the top N genes and use them for further evaluation.
Question: What are the differences between pcr primer and adapter? Are they all artificial sequence added in sequence amplification process? And can they both exist in raw sequence reads?

Thanks

Answer: Hello,

Both PCR primers and sequence adaptors are artificial DNA oligonucleotides, generally of known sequence. The first are used to prime DNA replication reactions. Those used to amplify specific DNA sequences *in vitro* are called PCR primers. Sequence adaptors are any kind of short DNA sequence serving the scope of fishing a (generally unknown) DNA sequence of interest for various purposes; they are used in a variety of techniques, and sometimes can take part in a DNA replication step (e.g. in a 5'RACE). Depending on experimental protocols and following bioinformatics, both PCR primers and sequence adaptors can end up in raw sequence reads (for instance, this may occur if you are sequencing a DNA library).
Question: Hi everyone,

I'm having a bit of trouble finding out how to use fimo to discover TFBS. The questions I had were:

How do you decide on good thresholds for PWM matching (--thresh)

How much do the backgrounds matter, should these be calculated from the sequences you're submitting?

Does how does Fimo handle masked fasta files? Is it better to submit hard of soft masked files?

Any other tips on good usage for this program would be much appreciated as I'm finding the documentation quite vague.

Answer: My response is probably not that great but since nobody else is answering this question:

you probably want to set a q-value of something like <0.05 or <0.01 since those cutoffs are pretty standard

backgrounds matter a lot, they should be calculated from the sequence that you are looking for motifs in.

how are you masking the files? with upper and lower cases or with X/Ns? you should be submitting sequence that you want to look for motifs in. If you are uninterested in motifs in repetitive regions, then you should mask them out (with Ns).

Lastly, have you tried using HOMER? The documentation on it is quite extensive and you might find it easier to use
Question: Hi,

I need help with awk commands, I have 4 samples in my vcf files, so field $10,$11,$12,$13 are the fields which have the genotype for each row, now I want remove the rows where in any of the any rows at least one sample is showing the genotype `./.` and want to print the rest in another vcf file, can this be done? Am not so familiar with awk substr. Any assistance? Below is the example of my vcf file, it does not have any header.

```
chr3 75787186 rs150410646 C T 53.89 . AC=4;AF=0.500;AN=8;BaseQRankSum=-4.341;DB;DP=424;Dels=0.00;FS=0.000;HaplotypeScore=2.2684;MLEAC=4;MLEAF=0.500;MQ=6.41;MQ0=371;MQRankSum=-3.553;QD=0.13;ReadPosRankSum=-1.007 GT:AD:DP:GQ:PL 0/1:63,21:80:48:48,0,127 0/1:25,5:29:21:21,0,64 0/1:142,41:174:10:10,0,94 0/1:95,31:120:6:6,0,120
chr3 75787576 rs141348932 A G 61.87 . AC=2;AF=1.00;AN=2;DB;DP=195;Dels=0.00;FS=0.000;HaplotypeScore=0.0000;MLEAC=2;MLEAF=1.00;MQ=4.17;MQ0=189;QD=0.69 GT:AD:DP:GQ:PL ./. ./. 1/1:68,22:86:9:87,9,0 ./.
chr3 75787583 rs144348996 A G 100.62 . AC=2;AF=1.00;AN=2;DB;DP=203;Dels=0.00;FS=0.000;HaplotypeScore=0.0000;MLEAC=2;MLEAF=1.00;MQ=4.33;MQ0=197;QD=1.12 GT:AD:DP:GQ:PL ./. ./. 1/1:65,25:86:12:126,12,0 ./.
chr3 75787584 rs151027881 C A 93.62 . AC=2;AF=1.00;AN=2;DB;DP=203;Dels=0.00;FS=0.000;HaplotypeScore=0.0000;MLEAC=2;MLEAF=1.00;MQ=4.33;MQ0=197;QD=1.04 GT:AD:DP:GQ:PL ./. ./. 1/1:64,26:86:12:119,12,0 ./.
chr3 75787620 rs145606249 T C 153.42 . AC=2;AF=1.00;AN=2;DB;DP=224;Dels=0.00;FS=0.000;HaplotypeScore=0.0000;MLEAC=2;MLEAF=1.00;MQ=4.38;MQ0=217;QD=1.70 GT:AD:DP:GQ:PL ./. ./. 1/1:52,38:86:18:179,18,0 ./.
chr3 75787728 rs111389701 C T 643.34 . AC=8;AF=1.00;AN=8;DB;DP=186;Dels=0.00;FS=0.000;HaplotypeScore=0.0000;MLEAC=8;MLEAF=1.00;MQ=10.21;MQ0=140;QD=3.46 GT:AD:DP:GQ:PL 1/1:0,32:32:3:28,3,0 1/1:0,23:23:9:82,9,0 1/1:0,82:82:51:503,51,0 1/1:0,49:49:6:55,6,0
```

I want to remove the rows where if any of the column $10,$11,$12,$13 is having `./.` no genotype then I want to eliminate those rows. Sorry for the formatting, I am not being able to get the correct format. Any suggestions?

Answer: You might take a look at [vawk][1], which is a tool from Aaron Quinlan's group and acts as an intelligent wrapper for awk on VCF.

[1]: https://github.com/cc2qe/vawk
Question: Hi,

I am trying to get the similar plot below for my data. I have VCF file from GBS analysis.
How can I draw this plot from VCF file.
Thanks.

![enter image description here][1]

[1]: /media/images/d27f72a9-9a0b-483d-8661-9de426c1

Answer: A similar plot can be generated with [CMplot][1] if you have data in such format.

[1]: https://github.com/YinLiLin/CMplot
Question: Hello

I am a beginner in Bioinformatics and I have a problem when I read some papers.

the paper name is SVDetect: a tool to identify genomic structural variations from paired-end and mate-pair sequence data and its context is related about a structural variants.

Anyway I can't figure out the meaning of strand orientation.

The word strand orientation comes from this sentence which is

> If one of the two ends of remaining pairs has an unexpected strand orientation, the cluster is annotated as a potential inversion.

I've already know that strand orientations means just direction of DNA(5'->3' or 3'->5')

but it seems strange considering the concept of strand orientation i knew. It can't match I think..

Help the meaning of strand orientation connected with the sentence above

Thank you

Answer: You stated that "strand orientations means just direction of DNA" but that is not correct.

DNA is made of two strands (who form a double helix). They run in opposite directions. Some proteins are encoded on the positive strand, some others on the negative one.

Each nucleotide is paired to its complement on the other strand; therefore you only need to record one strand.

Check [Wikipedia][1] or some MOOC about introductory Biology.

[1]: http://en.wikipedia.org/wiki/Dna
Question: Hi.

If you want to remove reads mapped on a certain chr, how would you do that?

I know 'samtools' can extract reads on a chr. Extraction of each chr and combining chr data I want would solve the problem, but I have multiple bam files to run it.

So, if anyone of you know how to REMOVE a single chr (e.g. chr1 only), please tell me.

Thank you!

Answer: See this post: https://www.biostars.org/p/128967/
Question: How to create a `Venn Diagram` `from data frame` and get the `list of common Genes` expressed in each combination?

I have a data frame of DEG, How can I create a `Venn Diagram` `from data frame` and get the `list of common Genes`?

Note: data frame has NA for some genotypes for some Genes, If all the genotypes are NA for a particular gene, that raw should be ignored.

Gene	Genotype1	Genotype2	Genotype3	Genotype4
AT1G17400	NA	NA	NA	NA
AT1G09420	NA	0.000800188	0.000116452	0.004017191
AT1G50930	NA	NA	NA	NA
AT1G65960	NA	NA	NA	NA
AT1G09400	NA	NA	NA	NA
AT1G09415	NA	NA	NA	NA
AT1G74730	NA	NA	NA	NA
AT1G75100	0.001639398	0.001578892	6.92E-05	NA
AT1G75100	0.001639398	0.001578892	6.92E-05	NA
AT1G75240	NA	5.60E-05	0.000235329	0.000162115
AT1G14920	NA	NA	NA	NA
AT1G14920	NA	NA	NA	NA
AT1G65510	NA	NA	NA	NA
AT1G75250	NA	NA	NA	NA
AT1G54410	NA	0.000113869	1.25E-05	NA

Answer: Question is very vague. What do these numbers mean? What do you consider expressed? Is anything that is not NA considered expressed?

Take a look at the Venn function here:
https://www.rdocumentation.org/packages/gplots/versions/3.0.1.1/topics/venn

It gives you an example on how to create Venn diagrams.

I'll go off the assumption that everything non-NA is considered expressed. Say you have everything stored in a dataframe named `data`. You need to supply lists of non-NA genes that belong to each of your four groups: Genotype1, Genotype2, Genotype3, Genotype4:

```
require(gplots)
Genotype1 <- data[!is.na(data$Genotype1),"Gene"]
Genotype2 <- data[!is.na(data$Genotype2),"Gene"]
Genotype3 <- data[!is.na(data$Genotype3),"Gene"]
Genotype4 <- data[!is.na(data$Genotype4),"Gene"]
input <- list(Genotype1=Genotype1, Genotype2=Genotype2, Genotype3=Genotype3, Genotype4=Genotype4)
venn(input)
```
Question: Hi

I want to do a python task which involves extracting the sequence from a multi-fasta file using a file of IDs for Mycosphaerella graminicola genome. I do not have data(multi-fasta file) for this task. Therefore, I downloaded from ENA using the following link. https://www.ebi.ac.uk/ena/data/search?query=Mycosphaerella+graminicola. Now I have 9 separate fasta files containing Mycosphaerella chr 3,chr 4,chr 5, chr 6, chr 10, chr 12, chr 13, chr 15, chr 19,chr 21 whole genome short gun sequencing.

My question 1)how to make a single multi-fasta file using these separate 9 fasta files. 2) Is the multi-fasta file available on web straightway?

Any Help would be highly appreciated!
Thanks

Answer: May be this also works:

$ awk "1" chr*.fa > out.fa
Question: I've read and been told that edgeR must take in counts only. I noticed that the `RPK` values are being fed into the `TMM` normalization procedure. Is this a correct usage assuming all of the assumptions?

Should this be used for downstream DGE analysis?

Note: I am no expert with these methods but I just wanted to ask the community

https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-018-2246-7

# calculate RPK
rpk <- (x[,2:ncol(x)]/x[,1])
# remove length col in x
x <- x[,-1]
# for normalization purposes, no grouping of samples
group <- c(rep("A",ncol(x)))
#EdgeR
x.norm.edger <- DGEList(counts=x,group=group)
x.norm.edger <- calcNormFactors(x.norm.edger)
norm.counts.edger <- cpm(x.norm.edger)

#GeTMM
rpk.norm <- DGEList(counts=rpk,group=group)
rpk.norm <- calcNormFactors(rpk.norm)
norm.counts.rpk_edger <- cpm(rpk.norm)

# Source:
# https://static-content.springer.com/esm/art%3A10.1186%2Fs12859-018-2246-7/MediaObjects/12859_2018_2246_MOESM4_ESM.docx

Answer: Technically, RPK values do not violate assumptions of TMM.

TMM is just a technique that tries to find the non-DE portion of the expression distribution by very liberally trimming off outliers. It doesn't matter what kind of expression units you are using.

However, RPK values do violate assumptions for DE analysis. So you cannot use it for downstream DGE.
Question: I have a file contain multiple gene >sequences. I want to separate into a file with gene ID perl or python script

Answer: This is a python script for splitting FASTA file into an individual file.

from Bio import SeqIO
import argparse

parser = argparse.ArgumentParser(description="Split the fasta file into individual file with each gene seq")
parser.add_argument('-f', action='store', dest='fasta_file', help='Input fasta file')
result = parser.parse_args()

f_open = open(result.fasta_file, "rU")

for rec in SeqIO.parse(f_open, "fasta"):
id = rec.id
seq = rec.seq
id_file = open(id, "w")
id_file.write(">"+str(id)+"\n"+str(seq))
id_file.close()

f_open.close()

To run above code, (save the above code in code.py file)

python code.py -f fasta_file

Note: You need to install Biopython module SeqIO to run this code.
Question: Hi all,

I am using mentha database for analyzing PPI networks. I have find the three centrality measures Degree,Closeness and betweenness centrality. Then I have sorted all three centralities from low to high. And now I am considering top 200 nodes. My network consists 16,000 nodes and 1,73000 interactions. Then I tried to find the common nodes between all the three centralities top 200 nodes. But the intersection between these three centrality results zero. That means the results are showing that there is no such common node between these three centralities. Is this kind of result is possible..?? And there is also no common node between closeness and betweenness and 4 common nodes between degree and closeness from top 200 nodes. So this kind of result is possible..?? This is a very rare result.

Answer: Although centrality measures tend to have positive correlations with each other, these measures all give different information so you shouldn't expect them to be high/low for the same proteins. The degree is the number of interactions that the protein makes, the closeness centrality measures how many "hops" in the network the protein is to other proteins and the betweenness centrality measures how often a protein is seen on the shortest path between other proteins. Imagine that the network represents relationships between people. If you want more people to buy you a drink, you should maximize your degree centrality i.e. have lots of friends. However, if you want to make sure that you know about every gossip making the rounds then you should maximize your betweenness centrality because having many friends doesn't help if they're isolated. Now if you're a lobbyist, you want to be able to talk to a wide range of people so you should maximize your closeness centrality.

There are other measures of centrality too that report on other properties of the nodes. Which centrality measures are relevant depends on the question you're asking about your nodes/network. For example, proteins with high degree but low closeness may form an isolated complex. Or high closeness and low betweenness could point to possible interaction "contaminants".
Question: I'm a bioinformatics student and I've never worked with the raw SNP data that you can get from 23andMe. What do people do with data like this? I'm not really looking for an analysis, but I would like to learn more about what a bioinformatician can do with a raw data set like this.

Here is the metadata and the first few lines as an example:

# This file contains raw genotype data, including data that is not used in 23andMe reports.
# This data has undergone a general quality review however only a subset of markers have been
# individually validated for accuracy. As such, this data is suitable only for research,
# educational, and informational use and not for medical or other use.
#
# Below is a text version of your data. Fields are TAB-separated
# Each line corresponds to a single SNP. For each SNP, we provide its identifier
# (an rsid or an internal id), its location on the reference human genome, and the
# genotype call oriented with respect to the plus strand on the human reference sequence.
# We are using reference human assembly build 37 (also known as Annotation Release 104).
# Note that it is possible that data downloaded at different times may be different due to ongoing
# improvements in our ability to call genotypes. More information about these changes can be found at:
# https://www.23andme.com/you/download/revisions/
#
# More information on reference human assembly build 37 (aka Annotation Release 104):
# http://www.ncbi.nlm.nih.gov/mapview/map_search.cgi?taxid=9606
#

rsid	chromosome	position	genotype
rs4477212	1	82154	AA
rs3094315	1	752566	AA
rs3131972	1	752721	GG
rs12124819	1	776546	AG

Answer: You can check the meaning of each SNP, given with the rs code. There are plenty snp databases where you can detailly read about the effects of each snp.
Question: I used Salmon to quantify RNA at the transcript level. Is there a way to map this data into a format that can be viewed in IGV?

A related question, i also tried another method of quantifying where i aligned the read with STAR and counted with GenomicAlignments. The counts from this method is easily double that of the count from Salmon. With some transcripts counted almost 8x to 10x more. I did a Pearson's correlation and it return as value of 0.95, so i'm guessing it would not significantly affect me if i am doing genome wide analysis. But i am still curious as to how this discrepancy came about.

Any input will be greatly appreciated

Regards

Answer: Hi,

I don't know about how to convert the data from salmon. I read, that you can generate pseudobam files with kalisto and maybe there is an tool for the output of salmon.

Concerning the STAR coverages, you get more reads, because it also counts gapped reads, which are not part of an already annotated transcript. Therefore, with salmon, where you only mapp the reads against annotated transcripts of a gene, you get a lot more reads, which also include splicing errors or not yet described alternative spicing events, where a gap can not be explained by already annotated transcript sequences. Salmon takes gapped reads into acoount, but only those who can be mapped to known transcripts. With gapped reads I mean gapped after alignment to the genome (like after STAR alignment). Thats why salmon an co are called pseudoalingers. They are fast and sufficient for well-annotated organisms. But you can't for example find new isoforms of a transcript.
Question: Hi all,

I have a question regarding coding regions in bacterial genomes: **can genes be encoded on both strands like this ?**

![enter image description here][1]

[1]: /media/images/fceb469c-9a40-47f9-969b-fdce9ddb

**If so, which one is considered as the sense strand ?**

Thank you so much for your help!

Answer: Yes genes can be encoded on both strands, not only in bacteria but also in viruses and eukaryotes.

The "sense" strand is relative to each gene. So it will be the top strand for the genes on the top strand (also called the Whatson or "+" strand) and vice-versa for the genes on the bottom strand (also called the Crick or "-" strand).