cannin/biostars_qa huggingface dataset:
Question: Hi all, I have just a very simple question about subsetting data.

I have a data frame with many columns, two columns are indicating the fold change after RNA-seq data analysis.
This two columns contain positive and negative values. I would like to keep the opposite values (negative in one column and positive in the other and vice versa) while removing those which are the same sign (negative and negative in the two columns or positive and positive).

Example:

df$1 df$2
1. 2
-1. 1
-1. 4
3. -2
6. -4
-6. -6
6. 1
-3. 7
......

Resulting dataframe :

df$1. df$2

-1. 1
-1. 4
3. -2
6. -4
-3. 7
......

Any help, highly appreciated, thanks!

Answer: Check the columns signs then filter, try this example:

# reproducible example data
df1 <- read.table(text = "col1 col2
1. 2
-1. 1
-1. 4
3. -2
6. -4
-6. -6
6. 1
-3. 7", header = TRUE)

# check if columns have the same sign, then filter
df1[ !(sign(df1$col1) == sign(df1$col2)), ]

# col1 col2
# 2 -1 1
# 3 -1 4
# 4 3 -2
# 5 6 -4
# 8 -3 7
Question: Hi,

I've recently started using Ensembl for help in designing a gene panel for NGS.

Each time I select a gene on Ensembl it comes up with different available transcripts for that gene.

What defines a transcript? Is it literally just a different version of the gene in different individuals (but if so, how come some of the transcripts appear to be so different in terms of bp length and amino acid number)

Any help in explaining would be most gratefully received!

Answer: A single gene can produce multiple different RNAs (i.e., transcripts). The actual transcript observed will depend on the tissue, developmental time point, and environmental/hormonal/etc. factors. Typically there's a single major transcript expressed in a given cell at a given time, but not always. For further details check out wikipedia.
Question: I have ChIP-seq data with peaks called by MACS in peaks.bed and summits.bed files. From here I would like to visualise the peaks on either the IGV or UCSC' Genome Browser. I understand that the bed files give me the positions and size of the peaks called by MACS, but how do I actually visualise the height / distribution of the peaks? Do I have to get the genome coverage for each base from my bam file, then filter for peak regions listed in the bed files, and then convert to a bedgraph, wig, or bigwig format?

Answer: 1\. Use `-w` and `-S` parameters, which gives you two 'wig' files: One for Control and another for TREATMENT. Wiggle (wig) file gives you the distribution and coverage (height) of peaks which you can upload on UCSC or IGV . Read more about the other parameters [here][1].

---Example:

macs14 -c -t -n OUTPUT_NAME -f BAM -g hs -w -S

2\. If you want to just show the coverage over predicted peak regions by MACS, use BEDtools to extract only the portion of 'wig' file overlapping with peak regions (BED file).

[1]: https://github.com/taoliu/MACS/blob/macs_v1/README.rst
Question: hello friends,

i have a list of GO terms, for example GO:0005975, how i can fine that which of these GO terms is related with cirus fruit senescence please???

Answer: For such a general question I would start with http://geneontology.org/ where you can begin a search based on "fruit senescence" (is there such a thing as cirus fruit?, did you mean citrus fruit?), as fruit senescence is more general and you're more likely to find results than requiring a particular type of fruit - unless deep knowledge already exists in that domain. However, this site will simply leave you over at http://amigo.geneontology.org/ where you can explore the relationships among terms and everything else. If you want a programmatic way of doing things - do what Deepak said.
Question: One can use CRAMtools to convert(compress) bam into cram. CRAMtools has options to allow for lossy and non-lossy compression (--lossless-quality-score and --lossy-quality-score-spec).

Samtools is also capable of compressing bam into CRAM but with limited parameteres, one just specifies the output to be cram (using option -C).
So my question is, do you know if the samtools' compression of bam into cram is lossy or lossless?

Thx

Answer: There is some sketchy documentation in the samtools man page (eg here: https://github.com/samtools/samtools/blob/develop/samtools.1#L1879 and further down for lossy names). These document options specified using --input-fmt-option or --output-fmt-option. Note it's not always obvious which is which without reading the man page. For example the required-field parameter is designed to avoid decoding certain components of a CRAM file, so when decoding CRAM or converting from CRAM to CRAM you can elect to simply drop certain fields. This permits total loss of quality if you so wished, or dropping of all aux tags. It's rather course and doesn't allow a lot of selection in tag types though. As this is a CRAM decode option it only works when reading CRAM. If you skip a data type it fills it out with the SAM default (so "*", 0, etc). Conversely the lossy names field is an output option which is used to discard query names when the read-pair are both within the same cram slice.

If you want more control over loosing tags and per-site degredation of qualities then maybe try Crumble. https://github.com/jkbonfield/crumble That too sadly needs documenting, or indeed writing up as a paper, but it's a matter of finding the time! (As always.)
Question: Assuming I would like to discuss something interesting about the gene when an adenine was substituted by a guanine (A>G as recorded in [rs268](https://www.ncbi.nlm.nih.gov/snp/rs268)), what is the appropriate nomenclature for me to do so? Is it appropriate to say **rs268 variant** when referring to the substituted sequence and **rs268 wildtype** when referring to the consensus sequence? What should I say if there are more than one substitutions recorded in the dbSNP entry?

Any input will be appreciated. Thanks in advance!

Answer: The correct terminology is "reference" and "alternative" (or "ref" and "alt"). The reference is whatever is in the reference genome, and is always listed first in allele strings. It is not necessarily the major, ancestral or wildtype allele, but it is in the reference genome so it is the reference allele. The alternative is whatever is not in the reference genome.
Question: I am new to Hi-C analysis.I am working with a test data set.can any one give me a clear work flow of Hi-C analysis and also can any one tell me about the maths behind the expected matrix and inter and intra chromosomal maps

Answer: HOMER contains a HI-C analysis tutorial:

http://biowhat.ucsd.edu/homer/interactions/index.html

http://biowhat.ucsd.edu/homer/interactions/HiCtips.html

A quote:

Hi-C is a beast of an assay. There is a lot of information in each
run, and the size and complexity of the experiments can stretch the
limits of your computer too.
Question: I'm using ballgown to analyze the data and got an error.

> pheno_data <- read.table(pheno_data_file, sep="\t", header=T)

> pheno_data

Samples	Group

Sample_male_1	Group1

Sample_male_2	Group1

Sample_male_3	Group1

Sample_female_1	Group2

Sample_female_2	Group2

Sample_female_3	Group2

> bg <- ballgown(dataDir = "ballgown", samplePattern="Sample", pData=pheno_data)

Wed Dec 21 10:36:10 2016

Wed Dec 21 10:36:10 2016: Reading linking tables

Wed Dec 21 10:36:15 2016: Reading intron data files

Wed Dec 21 10:36:32 2016: Merging intron data

Wed Dec 21 10:36:37 2016: Reading exon data files

Wed Dec 21 10:37:17 2016: Merging exon data

Wed Dec 21 10:37:24 2016: Reading transcript data files

Wed Dec 21 10:37:41 2016: Merging transcript data

Error in ballgown(dataDir = "ballgown", samplePattern = "Sample", pData = pheno_data) :
first column of pData does not match the names of the folders containing the ballgown data.

In addition: Warning message:

In ballgown(dataDir = "ballgown", samplePattern = "Sample", pData = pheno_data) :
Rows of pData did not seem to be in the same order as the columns of the expression data. Attempting to rearrange pData...

Execution halted

It seems that the order in phone_data should be exactly the same in the folder of ballgown.

Any help would be greatly appreciated..

Answer: hi
I recently got the same err.through reading the protocal,I find a solution:
To make sure that the file names match the
IDs in the phenotype file, run the following
R command:
> all(pheno_data$ids == list.
files("ballgown"))

I guess your pheno_data's order maybe not match the sample file order.
You can reorder your pheno_data_file's Samples as `list.
files("ballgown"))` to solve this problem.

Best
Question: Hi everyone! I am very new to the RNA-seq analysis and having quite a dilemma to design my experiment for the RNA-seq, in terms of the number of replicates required for each group of sample, with the consideration of budget constraint.

So here is the design of my experiment:
I am deriving different models of stem-like cells (or cancer stem cells, CSC) derived from HepG2 cell line, with a total of 3 models. So for each model, HepG2 cells were treated/cultured in different culture condition to enrich CSCs, meaning at the end of the day the total sample is 4 (1 parental HepG2 cells, 3 CSC-enriched models derived from HepG2 cells).

My questions here, if I do three replicates of these treatment/models, are they considered to be biological replicates? How are these considered biological replicate since they all come from the same cell line, just different treatment for the treated samples, aren't these technical replicates? Because I have read that technical replicates are unnecessary for RNA-seq but my concern is that whether the data will be reliable and publishable. But then again, is it right to do triplicates for RNA-seq from these samples? Plus, another concern are also considering the budgetary constraint of sequencing 12 samples, which are not affordable on our part.

Any thoughts? Sorry if my write-up is confusing, will really appreciate inputs from experts on this matter! Thanks!!!

Answer: You're correct that technical replicates are not necessary for RNA-seq, however biological replicates are always necessary for statistical analyses. Without replicates you cannot have any p-values.

A technical replicate in RNA-seq, in my experience, refers to taking one tube of RNA and then making multiple sequencing libraries from that tube. This is to control for errors in pipetting/PCR/anything technical during the library prep procedure. The kits available nowadays are so reproducible that this really isn't necessary. Biological replicates, however, are needed.

If I were designing the experiment, I would take the parental HepG2 cells, split them into three separate plates, then when the plates reach the appropriate confluence, I would treat them to generate the CSC-enriched models. Although there is an original single plate of cells, each of these biological replicates have undergone the CSC-enrichment separately, therefore any differences should be representative of this CSC-enrichment process (or biological differences).

You ***COULD*** get away with doing duplicates if sequencing 12 samples is absolutely unaffordable, some programs such as DESeq2 and edgeR are able to perform analyses with fewer replicates, but it really is not advisable. You can pool your 12 samples and run them down a single sequencing lane to save money, you'll have lower depth but you can always resequence and pool the raw fastq files if necessary.

I hope all of this is clear, please let me know if you have further questions.
Question: To identify the number of heterozygous variants in my .vcf file, I used the following linux command in vcftools:

$ vcftools --vcf SRR1611183.gatk.vcf --het

And returned the following:

INDV	O(HOM)	E(HOM)	N_SITES	F
NA12878	0	0.0	1101	0.00000

This is the same as the **'het_RA_count' => 1101** found in the `vcf-stats` feature; however, `vcf-stats` also lists **'het_AA_count' => 2**, which is not included in the `--het` command. If I'm seeking the total number of heterozygous variants should I just total them up (to 1103)? Perhaps het_AA is not truly heterozygous?

I'm open to a bcftools solution if its more straightforward.

Answer: You can try my 'independent' script and see what it returns: https://www.biostars.org/p/291147/#291167

That just counts the occurences of `0/1` `1/0` `0|1` and `1|0` per line. The total heterozyguos calls in the file will be the total of the final output column.

Another, more complex version of this is here: https://www.biostars.org/p/298361/#298464

Note that you should indeed be using BCFtools and not VCFtools. I believe that this advisory is even given by the VCFtools author (Adam).

Kevin
Question: Hello everyone,

I am fairly new to bioinformatics. I am currently learning the command line and was writing a very
simple script to align multiple files (Hepatitis C viral DNA sequencing) to reference genome using BWA, and then using samtools to align to the reference.

I seem to be getting a problem in the index step.

I keep getting the same error:

[E::hts_open_format] Failed to open file HCV242.sorted.bam
samtools index: failed to open "HCV242.sorted.bam": No such file or directory

Below is the whole script:

SAMPLE_LIST="HCV176 HCV242"
for SAMPLE_ID in $SAMPLE_LIST

do
#Index the aligner with the reference genome
bwa index -a is HCV1.fasta
# Align reads 1 of samples to ref genome
bwa aln HCV1.fasta $SAMPLE_ID.1.fastq >$SAMPLE_ID.1.sai
# Align reads 2 of samples to ref genome
bwa aln HCV1.fasta $SAMPLE_ID.2.fastq >$SAMPLE_ID.2.sai
# Merge and change into SAM format for both samples
bwa sampe HCV1.fasta $SAMPLE_ID.1.sai $SAMPLE_ID.2.sai $SAMPLE_ID.1.fastq $SAMPLE_ID.2.fastq > $SAMPLE_ID.sam
bwa sampe HCV1.fasta $SAMPLE_ID.1.sai $SAMPLE_ID.2.sai $SAMPLE_ID.1.fastq $SAMPLE_ID.2.fastq > $SAMPLE_ID.sam
# Change the SAM to BAM format for both samples
samtools view -bt HCV1.fasta $SAMPLE_ID.sam>$SAMPLE_ID.bam
# Sort the BAM files
samtools sort $SAMPLE_ID.bam -o $SAMPLE_ID.sorted
# Create the index file for the BAM files
samtools index -b $SAMPLE_ID.sorted.bam
done

I've tested the code and everything works until the index command. All the necessary files are generated as expected up until that step.
Any help much appreciated.

Answer: While very old versions of samtools sort would add the .bam at the end of the output file, I'm pretty sure yours does not. It should be easy enough to verify that you have $SAMPLE_ID.sorted, and not $SAMPLE_ID_sorted.bam
Question: I'm trying to figure out what the list of SNPs interrogated by [Peddy][1] is, and how one can manually specify the sites to look at, as stated in [their paper][2]:

> Given a VCF file and associated PED file describing the expected
> relationships and sex of the individuals in a sequencing study, peddy
> automatically conducts all of the tests described using the subset of
> 23,556 informative SNPs, thereby allowing the rapid detection of
> possible issues with individual samples. [...]
>
> We emphasize that the user may also specify their own selection of
> sites [...]

Does anyone have any experience with this?

Thank you.

[1]: https://github.com/brentp/peddy/blob/master/README.md
[2]: https://www.sciencedirect.com/science/article/pii/S0002929717300174?via%3Dihub

Answer: Here is the current list of sites: https://github.com/brentp/peddy/blob/master/peddy/1kg.sites

There is a grad student in the lab working on allowing any set of sites (so we can support hg38 and any genome build). I hope this will be out soon, but I think he has some work yet to do. You can track the progress here:
https://github.com/brentp/peddy/issues/33
Question: I have assembled a transcriptome using TRINITY with 100,000 transcripts. I have perform blastx against NR database (-outfmt 5). I want to remove redundancy in the assembled transcript before proceed to further processing.

1) How can I cluster them together into unigenes and remove redundancy?

2) Or how to select transcipts with longest sequences if they return same hit?

3) Any software or program suggestion for doing this?

4) Is it necessary for removing redundancy?

Thank you very much.

Answer: You will get transcripts after assembling short reads. Depending on alternative splicing,one gene can have more than one transcript.

Clustering or grouping puts all these similar sequences together and help to make a set of transcript for one gene.

On the other hand assembly process can create some transcripts that are not real (sequences/transcripts with more than 95% identity to another sequence) and clustering helps identifying them. Therefore, the flow goes like:

Assembly -> Transcripts -> clustering transcripts -> unigenes -> further downstream processing.

You can use [CD-HIT][1] for this purpose.

[1]: http://weizhongli-lab.org/cd-hit/
Question: Hello Forum,

What pipelines are available for exome sequencing analysis? Thanks for information.

Answer: In addition to Ram's excellent answer, have a look at [bcbio-nextgen][1] and [Brad Chapman's blog][2]. The former demonstrates some effective methods for alignment and processing, the later describes a good number of comparisons that you should be familiar with if you want to really get into variant calling pipelines.

[1]: https://bcbio-nextgen.readthedocs.org/en/latest/
[2]: http://bcb.io/
Question: Hello,

Three sources suggest different read lengths -- how can I definitively find the read length?

1. A colleague suggested that this bam file has reads 40 base pairs long,
2. samtools view seems to suggest the reads are 44 base pairs long, and
3. [Tablet][1] seems to suggest most reads are 40 base pairs long but shows a couple less than 40 (e.g. 36, 38, and 39) base pairs long.

For samtools, I'm running the following command, based off of https://www.biostars.org/p/65216/:

```
$ samtools view foo.bam | head -n 1000000 | cut -f 10 | perl -ne 'chomp;print length($_) . "\n"' | sort | uniq -c
1000000 44
$ samtools view foo.bam | head -n 1
HWI-ST673_0087:1:1103:11279:92479#NNNNNN/1 0 chrM 1 255 4M4I36M * 0 0 GATGGATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCT bbbeeaceggggfeghiiiiiiiiiiihiiihiiiiihfhhhhi AS:i:-23 XN:i:0 XM:i:1 XO:i:1 XG:i:4 NM:i:5 MD:Z:3C36 YT:Z:UU
```

With tablet, I see the following:

I assume that I'm just misinterpreting the output. I don't have any formal background in bioinformatics (I am a software developer). If anyone could help set me straight here, I'd appreciate it, or maybe suggest some good relevant reading. Do bam files always have uniform read lengths?

Thanks for the help!

#### Update: was looking at different file in tablet -- discrepancy resolved!

[1]: http://ics.hutton.ac.uk/tablet/

Answer: The sequence reported in the bam file is what I would trust for the read length. This may be different from the alignment length, due to clipping, insertions, or deletions (see [here][1] for more). The visualizers you use may only plot the aligned portion of the read, which could lead to the differences that you're seeing.

There are no rules about read lengths in bam files, since you can put anything into them and they may represent combinations of various experiments. My experience is that a single Illumina sequencing run will have constant read length, but if your files represent multiple runs or have preprocessed data (maybe via adapter or poly-A trimming), that assumption won't hold.

[1]: http://davetang.org/wiki/tiki-index.php?page=SAM#Clipped_alignment
Question: I need to calculate the median in the fastq file, how can I do this with awk or bioawk?

Answer: with seqkit

$ seqkit stats -a test.fq

In output, Q2 is median value. If you want to extract fastq name, and median, try

$ seqkit stats test.fq -aT | cut -f1,10

with awk and st (available in brew, conda):

$ awk 'NR % 4==2 {print length ($1)}' hcc1395_tumor_rep1_r1.fastq |st --summary

min	q1	median	q3	max
151	151	151	151	151

$ awk 'NR % 4==2 {print length ($1)}' hcc1395_tumor_rep1_r1.fastq |st --median
151

Ministat is available in system repos for many GNU-Linux and BSD distros and you can try:

$ awk 'NR % 4==2 {print length ($1)}' hcc1395_tumor_rep1_r1.fastq | ministat -A |tr -s " "| cut -f5 -d " "

Median
151
Question: Hello Biostars,

What are an approximate quantity of human genes, transcripts (coding transcripts), proteins and metabolites? What about a model plant like Arabidopsis?

Thanks

Answer: # Human

- [GENCODE ](http://www.gencodegenes.org/stats/current.html)

- [HGNC](http://www.genenames.org/cgi-bin/statistics)
Question: I have a set of differentially expressed genes after analysis of RNAseq data. How can I sort them into protein coding genes and miRNA, ncRNA?
I have used HISAT2-STRINGTIE-BALLGOWN for analysis.
Reference genome: GRCh38 with annotation from Ensembl release 84.

Answer: Colsult our ancestors wisdom:

[Question: How One Can Separate Noncoding From Coding Ucsc And Ensembl Transcripts][1]

[Question: Identify lncRNA in list of Ensembl ID's][2]

[Question: How to filter out the non-coding genes?][3]

[Question: How to separate htseq-count table to coding RNAs and non-coding RNAs][4]

[Question: BiomaRt: protein coding transcript identification][5]

[1]: https://www.biostars.org/p/72009/
[2]: https://www.biostars.org/p/199646/
[3]: https://www.biostars.org/p/101643/
[4]: https://www.biostars.org/p/159569/
[5]: https://support.bioconductor.org/p/62441/
Question: Are there any public datasets that have Unique molecular identifiers of some sort in the reads? It would be interesting to see or test that data, but I am not familiar with any available data.

Also, I tried the data in this paper https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4907374/. The data sets seem to be already trimmed of UMIs, as far as I can tell, or I might be missing something.

Does anyone have any advice on datasets?

Thanks in advance!

Answer: Try the datasets referenced in this paper:

https://www.ncbi.nlm.nih.gov/pubmed/28100584
Question: I'm trying to get a regex to work with rename; I've tried the approach of similar answered questions here but couldn't get the results I wanted.

The files are named as such:

SR1_S90_L001_R1_001.fastq.gz
SR1_S90_L001_R2_001.fastq.gz
Rinc_S96_L001_R1_001.fastq.gz
Rinc_S96_L001_R2_001.fastq.gz

And I would like to retain only the information prior to the first underscore and the _R1_ or _R2_ tags, like this:

SR1_R1_.fastq.gz
SR1_R2_.fastq.gz
Rinc_R1_.fastq.gz
Rinc_R2_.fastq.gz

Thanks in advance!

Answer: `rename -n 's/(\w_).*_(R[0-9])_.*(.fastq.gz)/$1$2$3/' *.fastq.gz` or `rename -n 's/(\w+_)\w+_\w+_(\w._)\w+(.\w+)/$1$2$3/' *.fastq.gz`

-n runs the command in dummy mode and it is distro specific. Check the available for options for rename on your distro. -n option is available on ubuntu 18.04 and remove -n for final conversion.