cannin/biostars_qa huggingface dataset:
Question: I have a data frame with the differentially expressed genes from EdgeR, Now I am trying to make a volcano plot of it but I want to see only selected genes that are of interest to me to be labelled on the volcano plot.
My data frame looks like this

head(results)

Gene Fold pvalue FDR sig
1 ADORA2A 10.273854 1.164636e-28 2.234471e-24 FDR<0.05
2 IL23A 8.132293 4.554177e-28 4.368822e-24 FDR<0.05
3 IL1A 8.430078 6.768343e-27 4.328581e-23 FDR<0.05
4 CXCL6 7.102900 3.299464e-23 1.582588e-19 FDR<0.05
5 CCR7 8.950486 9.111421e-23 3.496235e-19 FDR<0.05
6 IL18R1 6.759646 7.283440e-22 2.329001e-18 FDR<0.05

I tried the below code to generate the volcano plot, which is generated successfully and the selected genes are marked too but not with their exact name but with some random number.

results$genelabels <- ""
results$genelabels <- ifelse(results$Gene == "IL23A"
| results$Gene == "IL1A"
|results$Gene == "IL6"
|results$Gene == "CD80"
|results$Gene == "CD86"
|results$Gene == "NFKB"
|results$Gene == "BAFT2", TRUE,FALSE)
ggplot(results) + geom_point(aes(Fold, -log10(FDR),col=sig))+ geom_text_repel(aes(Fold, -log10(FDR)),label = ifelse(results$genelabels == TRUE, results$Gene,""), box.padding = unit(0.45, "lines"),hjust=1) + theme(legend.title=element_blank(),text = element_text(size=20))+ scale_color_manual(values=c("red", "black"))

Can somebody help me out with what is wrong with the code that is does not show the exact name of the gene on volcano plot?

Answer: Hi [saamar.rajput][1],

See the answer from [benformatics][2] (credits goes to him), the key is to make sure your "Gene" is in character type instead of factor.

```
> results <- read.delim("example.txt", header = TRUE, stringsAsFactors = FALSE)
> head(results)
Gene Fold pvalue FDR sig
1 ADORA2A 10.273854 1.164636e-28 2.234471e-24 FDR<0.05
2 IL23A 8.132293 4.554177e-28 4.368822e-24 FDR<0.05
3 IL1A 8.430078 6.768343e-27 4.328581e-23 FDR<0.05
4 CXCL6 7.102900 3.299464e-23 1.582588e-19 FDR<0.05
5 CCR7 8.950486 9.111421e-23 3.496235e-19 FDR<0.05
6 IL18R1 6.759646 7.283440e-22 2.329001e-18 FDR<0.05
> str(results)
'data.frame':	6 obs. of 5 variables:
$ Gene : chr "ADORA2A" "IL23A" "IL1A" "CXCL6" ...
$ Fold : num 10.27 8.13 8.43 7.1 8.95 ...
$ pvalue: num 1.16e-28 4.55e-28 6.77e-27 3.30e-23 9.11e-23 ...
$ FDR : num 2.23e-24 4.37e-24 4.33e-23 1.58e-19 3.50e-19 ...
$ sig : chr "FDR<0.05" "FDR<0.05" "FDR<0.05" "FDR<0.05" ...

> results$genelabels <- ""
> results$genelabels <- ifelse(results$Gene == "IL23A"
+ | results$Gene == "IL1A"
+ | results$Gene == "IL6"
+ | results$Gene == "CD80"
+ | results$Gene == "CD86"
+ | results$Gene == "NFKB"
+ | results$Gene == "BAFT2", TRUE, FALSE)

> ggplot(results) +
+ geom_point(aes(Fold, -log10(FDR), col = sig)) +
+ geom_text_repel(
+ aes(Fold, -log10(FDR)),
+ label = ifelse(results$genelabels, results$Gene, ""),
+ box.padding = unit(0.45, "lines"),
+ hjust = 1
+ ) +
+ theme(legend.title = element_blank(), text = element_text(size = 20)) +
+ scale_color_manual(values = c("red", "black"))
```

![volcano][3]

[1]: https://www.biostars.org/u/20971/
[2]: https://www.biostars.org/p/381048/
[3]: https://i.ibb.co/X5xVPth/Rplot.png
Question: This is a basic question. I want to rearrange a table like this:

```
A 1
B 2
B 3
B 4
C 5
C 6
```

to:

```
A 1
B 2,3,4
C 5,6
```

Could anyone gives a hint. Either with R or python?

Answer: How about with bedtools:

$ cat t.txt
A 1
B 2
B 3
B 4
C 5
C 6

$ bedtools groupby -i t.txt -g 1 -c 2 -o collapse
A 1
B 2,3,4
C 5,6
Question: Hello all,

I am a newbie in bioinformatics and would like to process a Trio (Parents and Index) in order to understand the concept. Here is my question: How can I analyse a WGS Trio in order to conduct paternity testing? Or is it possible to do this?

The question may sound stupid but I am very curious about it.

Many thanks!

Answer: It is definitely *possible* to assess paternity from whole genome sequence (WGS) data. Paternity can probably be established with as little as a few dozen or maybe hundreds of well-chosen single nucleotide polymorphisms (SNPs). If you have decent WGS data you can expect to genotype millions of SNPs. So, paternity assessment would be very confident from such data. What data do you have available? Assuming that you have raw sequence data (e.g., fastq or unaligned bam files) you will first need to align to an appropriate reference genome.

There are several online tutorials to give you the general idea:

- https://icb.med.cornell.edu/wiki/index.php/Elementolab/BWA_tutorial
- https://www.biostars.org/p/41951/

Note. Both the above tutorials are a little out of date. Current best practice would be to use bwa mem (available with current bwa installations). See http://bio-bwa.sourceforge.net/bwa.shtml

Once you have aligned your data you will probably want to mark duplicate reads and perform base quality score recalibration (BQSR). For some sample commands taking you through bwa mem alignment, duplicate marking, and BQSR see here: http://pmbio.org/module%202/0002/01/31/Alignment/

Next, you will want to run GATK variant caller. For a trio analysis I suggest you try running GATK HaplotypeCaller in GVCF mode and then performing joint genotyping. See here for a tutorial on this topic:
https://gatkforums.broadinstitute.org/gatk/discussion/7869/howto-discover-variants-with-gatk-a-gatk-workshop-tutorial

This is all explained in great detail in the excellent GATK Best Practices for Variant Discovery workshops organized by the Broad. See
https://drive.google.com/drive/folders/1U6Zm_tYn_3yeEgrD1bdxye4SXf5OseIt

Finally, assuming you get through the above. You should have a VCF with genotype calls for millions of SNPs for your trio. You then need to look at SNP genotype concordance between individuals in your trio to estimate kinship. This is itself a complicated area of research that I am not very familiar with. But, paternity should be one of the simpler relationships to prove. I believe the KING tool is popular for this and could take the above VCF as a starting point.

http://people.virginia.edu/~wc9c/KING/manual.html
Question: As a bioinformatics team grows there is huge benefit from sharing code, particularly with strong coding guidelines that promote read-ability and reuse.

But can we really pick common technology, infrastructures shared with our colleagues? Generic scripting languages like PERL, Python and Groovy are great, and each come with different strengths for common file and data marshalling or common mathematically/statically data field manipulations.

Other less generic tools, like R, KNIME, Galaxy take us down deeper applications for data analysis without becoming science application specific. And often it is specific tools for the instrumentation that is more efficient for targeted work. What is a bioinformatics team to do?

Thoughts?

Answer: Use the best tool for the job? Have people around that know Perl or Python and R, are very efficient under Linux and can learn to use new tools fast. The problems that are found within the field of bioinformatics are very diversified and you won't come close to anything that could be regarded as 'the best generic scripting tool'. Cheers
Question: Hello guys,

I need any tool or packages for integrated analysis of miRNA and mRNA expression data? More specific I have time course expression data from microarray experiment.

Thanks in advance

Answer: It depends on your definition of integrated analysis. The most simple (and also least powerful) way of integration is to do differential analysis on mRNA and miRNA between your timepoints seperately and afterwards compare the differential results. You can look if annotated mRNA-miRNA pairs according to the miranda database change over time in opposite directions. However, a much more powerful approach is to use regression of miRNA on mRNA expression estimates to look for significant negative interactions between any miRNA-mRNA pair. You don't need a special package, linear regression is in the base functionality of R. The advantage of the latter approach is explained by the following example. Supposed you have a set of 50 control and 60 tumor samples. In all samples, the miRNA expression and mRNA expression, which are known to be interacting, are similar except for 1 tumor. This tumor has low expression of the mRNA and high expression of the miRNA. So you see silencing of the gene that is most likely caused by the increased expression of the miRNA. With regression analysis you will find this very easily. However if you did differential expression analysis between tumor and control, both miRNA and mRNA would not have been significant because the one tumor sample does not make the means of the tumor and control cohorts different.
Question: Hello,
I'm analyzing some ChIP-seq data (bam files) through deepTools' bamCoverage and bamCompare. I'm new to biological data analysis and I don't understand the difference between scaling (option --scaleFactorsMethod in bamCompare) and normalizing (option --normalizeUsing).

BamCompare doesn't allow the usage of both (scaling the 2 compared bam files and individually normalizing the bam files).

My questions are:

- What do they mean with scaling? When should I use it?

- What do they mean with normalization? When should I use it?

Which one accounts for sequencing depth?

Thanks,

Lorenzo

Answer: Normalizing put values on a (somewhat) universally coherent value range. For example, `--normalizeUsing RPGC` will normalize the sample(s) to 1x coverage before proceeding. Scaling is done when you have 2 samples and needs to be done to make them comparable (i.e., it corrects for differences in sequencing depth).

- Using normalization whenever you need to compare multiple samples, or when it's the relative change in signal that matters (e.g., ChIP-seq)
- Scaling alone should be used when you need to compare samples but need the output to mimic the raw read densities. This is often useful for things like ribosomal-profiling.

For ChIP-seq stick to normalization.
Question: I want to extract some sequences using ID from a multifasta file.

Using perl is not possible because it gave an error when indexing the database. Maybe because of it's size? Is there any way to this differently? thanks for your help!

Answer: How about samtools faidx? It seems to be less well-known, but I think it should work for you.

If you have these reads in test.fa:

>one
AAAAA
>two
CCCCC
>three
GGGGG

Then just run:

samtools faidx test.fa

This will create an index file that will allow you to do sublinear queries for each name lookup (faster than the awk/grep solutions above and easier than rolling your own script).

For example:

samtools faidx test.fa two

Will output:

>two
CCCCC

It works for multiple regions (samtools faidx test.fa one two three) and that lets you play tricks like:

xargs samtools faidx test.fa < names.txt
Question: After running cufflink for two samples of RNA-Seq data, I used cuffmerge and executed the following command

cuffmerge -p 8 -g / -s / /assemblies.txt

Where `assemblies.txt` contains the `transcripts.gtf`

But I am getting following error

```
[Sun Jul 13 14:44:46 2014] Beginning transcriptome assembly merge
-------------------------------------------

[Sun Jul 13 14:44:46 2014] Preparing output location ./merged_asm/
[Sun Jul 13 14:46:04 2014] Converting GTF files to SAM
[14:46:04] Loading reference annotation.
[14:47:10] Loading reference annotation.
[14:48:06] Loading reference annotation.
[14:49:23] Loading reference annotation.
[14:50:29] Loading reference annotation.
[14:51:44] Loading reference annotation.
[Sun Jul 13 14:52:46 2014] Quantitating transcripts
Warning: Your version of Cufflinks is not up-to-date. It is recommended that you upgrade to Cufflinks v2.2.1 to benefit from the most recent features and bug fixes (http://cufflinks.cbcb.umd.edu).
Command line:
cufflinks -o ./merged_asm/ -F 0.05 -g / -q --overhang-tolerance 200 --library-type=transfrags -A 0.0 --min-frags-per-transfrag 0 --no-5-extend -p 8 ./merged_asm/tmp/mergeSam_file3DKFYb
[bam_header_read] EOF marker is absent. The input is probably truncated.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
File ./merged_asm/tmp/mergeSam_file3DKFYb doesn't appear to be a valid BAM file, trying SAM...
[14:52:49] Loading reference annotation.
Error parsing strand (?) from GFF line:
IWGSC_CSS_1DS_scaff_731014 . repeat_region 1 174 . ? .Name=trf;class=trf;repeat_consensus=ATTGGTATAGAACGCATGAAGAAACTCCATACAGATGGATCTTTAGACTCACTCAATTATGAAAAAATTGAGACATGCAAACCATGTCT;type=Tandem repeats
[FAILED]
Error: could not execute cufflinks
```

Answer: You are using the wrong input file. You downloaded the GFF3 instead of GTF file. The correct file is: ftp://ftp.ensemblgenomes.org/pub/plants/release-22/gtf/triticum_aestivum/Triticum_aestivum.IWGSP1.22.gtf.gz

This is your GFF3 file:

```
##gff-version 3
IWGSC_CSS_5BS_scaff_1034127 . repeat_region 61 200 . + . Name=gnl|TREP|TREP3026;class=Unknown;repeat_consensus=N;type=Unknown
IWGSC_CSS_3AL_scaff_747250 . repeat_region 2 200 . + . Name=gnl|TREP|TREP765;class=Unknown;repeat_consensus=N;type=Unknown
IWGSC_CSS_3B_scaff_7107049 . repeat_region 1 200 . + . Name=gnl|TREP|TREP232;class=Unknown;repeat_consensus=N;type=Unknown
```

This is the GTF file:

```
IWGSC_CSS_6DL_scaff_127793 protein_coding exon 526 645 . + . gene_id "Traes_6DL_7FFFE462C"; transcript_id "Traes_6DL_7FFFE462C.2"; exon_number "1"; seqedit "false";
IWGSC_CSS_6DL_scaff_127793 protein_coding CDS 574 645 . + 0 gene_id "Traes_6DL_7FFFE462C"; transcript_id "Traes_6DL_7FFFE462C.2"; exon_number "1"; protein_id "Traes_6DL_7FFFE462C.2";
```
Question: The term HSP has been used in many research papers while describing BLAST output. What are High Scoring Pairs (HSP)? I'm not very successful in Googling.

Answer: http://www.ncbi.nlm.nih.gov/books/NBK62051/

> A High-scoring Segment Pair (HSP) is a local alignment with no gaps
> that achieves one of the highest alignment scores in a given search.
Question: I have a bam file that I want to split into paired end fastq files per lane. Can somone please suggest me how it can be done?

Answer: Using the example data you have above it can be demultiplexed into lane specific files using `demuxbyname.sh` from BBMap suite as follows. (use `.gz` extension to keep files zipped). Input here is interleaved (as in your example). If someone reaches this solution in future and have their inputs in two separate files then specify them with `in1=R1` and `in2=R2`.

$ demuxbyname.sh -Xmx10g in=file.fq delimiter=: column=2 out=%_#.fq

Example output

$ more 5_1.fq
@HS2000-1015_160:5:2306:10070:71746/1
AGAAGACGAGCACAGTGAGGAAGGAGATGAGCAGGCAGGGGATGATGAGGTTGATGGTGTAGAACAAGGGCAGGCGCCGGATGTACAGCGAGTATGTGAT
+
CCCFFFFFHHHHHIICGGHHHAFGEFGIJJJIIJIIGIIGG;FGHEHHGI=CHCEHGHE?DFFDFECEDDDDDDDDDBDDBD@ACDAACDBD5@>
@HS2000-1015_160:5:2113:11446:94436/1
AATTGCTTGAACCCGGGAGGCAGAGGTTGCAGTGAGCCAAGATCGCACCACTGCAGTCCAGCCTGGGAAACCAGAGCAAAACTTGGTCTCAAAAAAAAAA
+
CCCFFFFFHHHHGHGFGFHJIBF1DC<@;)9D*099?F;8=(BFF5(;'.=CEHFBDC7.>;?C@BBB35>(AC?B>4+(+>A:(:@(2)5@B#

$ more 5_2.fq
@HS2000-1015_160:5:2306:10070:71746/2
GAACCTCAAGGACTATTGGGAGAGCGGCGAGTGGGCCATCATCAAAGCCCCAGGCTACAAACACGACATCAAGTACAACTGCTGCGAGGAGATCTACCCC
+
@CCFFFFDHGHFHIJJJJJJGGGIIJJIGHI@FHGIIGHHEFGHHFFFFFBCDDDDDCDDDDDDDD;@BDCCDACDD@>ACCDDDDBDB;?==AABC95<(5(5309@D########################
Question: Hello,

I need help using the ***CommandLine mode*** of MEGAN5.

I want to process a large set of Blast_comparisons using MEGAN5 under Linux. I use Megan5 routinely and am happy with it, but for my large number of samples I don't want to repeatedly click on the import options every time manually.
Instead, I want to use the commandline options of MEGAN to automize this.

However I am having great trouble figuring out exactly how to do this.
My data is in tabular BLAST format with GI numbers as sequence-IDs.
I do not need to import/embed the actual sequences in my MEGAN-file.
I use the mapping files provided by the MEGAN developers to map the GI numbers in my Blast-data to Taxon-information

I tried this set of commands in a scriptfile ('MEGAN_COMMANDS.txt') for simply importing my blast_data into MEGAN5:

load taxGIFile='/usr/local/share/db/MEGAN5_mapping/gi_taxid_nucl.bin'
set taxUseGIMap=true
import blastFile='/home/me/MYBLAST.tab.out' meganFile='/home/me/MYBLAST.tab.out.rma' maxMatches=100 minScore=1.0 maxExpected=0.01 topPercent=10.0 minSupport=1 minComplexity=0.0 useSeed=false useCOG=false useKegg=false usePfam=false paired=false useIdentityFilter=false textStoragePolicy=Embed blastFormat=BlastTAB mapping='Taxonomy:BUILT_IN=true,Taxonomy:GI_MAP=true,SEED:GI_MAP=false,KEGG:GI_MAP=false'

I tried supplying these commands to MEGAN5 like this:

MEGAN -g -v -E < MEGAN_COMMANDS.txt

And like this:

MEGAN -v -c MEGAN_COMMANDS.txt

But in both cases MEGAN just opens a blank screen and then just hangs (I have to abort it using the "kill" command in a new shell).
If I import the data by clicking manually it imports the data correctly in a matter of only a few seconds

What am I doing wrong here?

Answer: Finding [this thread][1] on SeqAnswers finally helped me find the problem!
In fact my problem was rather more simpler (and stupid) than in the mentioned thread:
Each command line in the command-file ("MEGAN_COMMANDS.txt" in my case) has to end with a semicolon. After appending semicolons to each line it worked!

Now I can finally automize my MEGAN analyses (GUI's may be nice sometimes, but just not practical for repetitive tasks)

[1]: http://seqanswers.com/forums/showthread.php?t=43535
Question: The GSEA manual states that the signal2noise ranking metric should be used only in case of 3+ samples per condition.

What is the reason for this restriction?

I do understand that the denominator has standard deviations for the input values and it may be more meaningful in case of 3+ samples but how is it unsuitable in case of 2 samples?

Answer: I guess it is partially arbitrary, but the confidence intervals improve a lot with 3 samples, compared to 2. Suppose a SD of 1, then you have the following [confidence intervals][1]:

N=2	95% CI of the SD	0.45 to 31.91
N=3	95% CI of the SD	0.52 to 6.28
N=4	95% CI of the SD	0.57 to 3.73
N=5	95% CI of the SD	0.60 to 2.87

[1]: https://www.graphpad.com/quickcalcs/CISD1.cfm
Question: I have WES/WGS data on a normal embryonic stem cell line and I want to cell any possible variants in it. I worked with the GATK germline variation calling and I got the results that I should've gotten if I had a germ cell.
Is that normal? Can you use germline variation calling on a normal embryonic stem cell line?

Answer: I think maybe you're confused about what is meant by "germline". When talking about "germline variant calling" we're not referring to the cell type, we're referring to where the variant came from. In this case "germline" means that the variant was inherited from a parent and should be found in every cell of the individual. This is opposite to "somatic" which means the mutation occurred somewhere in that individual's body and will be shared only by a small number of cells, often a tumour. The process of "germline variant calling" is, therefore, identifying all the variants in that individual that differ from the reference genome, and the process is identical whatever cell type you use for this.
Question: Hi,

I have downloaded Geuvadis genotypes from:

https://www.ebi.ac.uk/arrayexpress/files/E-GEUV-1/GEUVADIS.chr22.PH1PH2_465.IMPFRQFILT_BIALLELIC_PH.annotv2.genotypes.vcf.gz

And I am trying to recreate vcf files with only the European samples using bcftools, but when I try

bcftools view --samples-file 373_sampleIDs.tab GEUVADIS.chr22.PH1PH2_465.IMPFRQFILT_BIALLELIC_PH.annotv2.genotypes.vcf.gz

I get this error:

[W::vcf_parse_format] FORMAT 'PP' is not defined in the header, assuming Type=String
[W::vcf_parse_format] FORMAT 'BD' is not defined in the header, assuming Type=String
Undefined tags in the header, cannot proceed in the sample subset mode.

Trying to fix it with bcftools reheader does not work:

bcftools view -h file.vcf > header.txt; bcftools reheader -h header.txt file.vcf > fixed.vcf

I can find information on the FORMAT "PP", but not on the FORMAT "BD".

##FORMAT=

Would anybody be able to help?
Thanks!

Answer: the following command seems to fix the error. But I've no idea about the description of the new FORMAT. There are also some extra spaces in the INFO column...

wget -O - "https://www.ebi.ac.uk/arrayexpress/files/E-GEUV-1/GEUVADIS.chr22.PH1PH2_465.IMPFRQFILT_BIALLELIC_PH.annotv2.genotypes.vcf.gz" |\
gunzip -c |\
awk '/^#CHROM/ {printf("##FORMAT=\n##FORMAT=\n");} {print}' |\
sed 's/ damaging/_damaging/g' > fixed.vcf
Question: Hello All,

I am working on TCGA lung cancer data , I want to compare Average expression of a set of gene (my interested set of gene) in normal and tumor samples. I am wondering that the average expression of these gene in normal and tumor samples are very simialr in normalized log2 data , Fig1, (LUAD.uncv2.mRNAseq_RSEM_normalized_log2.txt), but it is different in normalized Z_score data, Fig2, (LUAD.uncv2.mRNAseq_RSEM_Z_Score.txt).

Fig 1, when using (LUAD.uncv2.mRNAseq_RSEM_normalized_log2.txt) data

![Fig 1][1]

Fig 2, when using (LUAD.uncv2.mRNAseq_RSEM_Z_Score.txt) data
![Fig2][2]

PS: x-axis, same order of genes

So, please help me, which input data should be approprite for this type of comparision.

Thank you...
[1]: http://s33.postimg.org/9pc7pa9q7/image.png
[2]: http://s33.postimg.org/3qccl35q7/image.png

Answer: Here, Z_Score means:

Z_Score = (expression in single tumor sample) - (mean expression in all tumor samples ) / (standard deviation of expression in all tumor samples)

That's why the Z-score for cancer group is very small in your Figure 2.

And I think for the calculation of Z score for normal samle is something like this way:

Z_Score = (expression in single normal sample) - (mean expression in all normal samples ) / (standard deviation of expression in all normal samples)

That means this curve only show the fluctuation of the gene expression in that group

I am pretty sure that you should use the data in Figure 1.
Question: I have a plot,where y-axis are obervations and x-axis is time. I have to find the maximun (x,y point) from the plot. I used max(), but I only figured out the obervation value but not the time point where this max obervation value is. Could someone help me please?

Answer: Try this example, using base plot:

# dummy data
df1 <- mtcars[1:10, 1:2]

# plot
plot(df1$mpg, df1$cyl)

# get the max on y, and matching x value
df1Max <- df1[ which(max(df1$cyl) == df1$cyl), ]

# mark the max point on the plot
points(df1Max$mpg, df1Max$cyl, col = "red", pch = 19)

# add vertical line
abline(v = df1Max$mpg, col = "blue", lty = "dashed")

![enter image description here][1]

[1]: https://i.imgur.com/jYTqdPm.jpg
Question: Hi,

I was wondering what happens if you have the libtype wrong for Salmon? Or why it requires that you supply it? I wasn't sure if the first read of my paired-end reads was forward (ISF) or reverse (ISR), so I ran it both ways, and the mapping efficiencies were identical between the runs.

Thanks,
Matt

Answer: Hi [Matthias](https://www.biostars.org/u/30561/),

You can find a description of the behavior you're seeing, and why it's generally expected [here](https://github.com/COMBINE-lab/salmon/issues/67). Also, I'll note that if you update to the latest release of Salmon (currently v0.7.0), available [here](https://github.com/COMBINE-lab/salmon/releases/tag/v0.7.0), you can use the `A` (automatic) library type to have Salmon automatically infer the most-likely library type for you.