cannin/biostars_qa huggingface dataset:
Question: Hi guys,
I've the following dataframe df:

IID PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9
1 11:00223.CEL 0.00229647 -0.000423608 0.001480000 1.02983e-03 0.00418171 -0.00550339 0.003826840 0.002836460 -0.001210240
2 9399.CEL 0.00213518 -0.000734612 0.000965396 9.84589e-04 -0.00135706 0.00396061 0.006356090 0.000639752 -0.000220536
3 shckanc.CEL 0.00225502 -0.000971542 0.001553290 2.13150e-03 -0.00277924 0.00291143 0.007620090 -0.002271640 -0.003364920
4 787323A1.CEL 0.00292725 -0.000399576 0.001340910 6.99230e-06 -0.00299630 -0.00769350 -0.007620050 0.002889230 0.005459370
5 90w3jc.CEL 0.00228869 -0.000201784 0.001291800 8.68436e-04 -0.00194812 -0.00554044 -0.006723900 0.001160370 0.005532150
6 olsj909.CEL 0.00224916 -0.000530798 0.000677401 -1.13087e-04 -0.00132783 -0.00814353 -0.000705561 0.000352494 0.000671816

Considering the first column IID, I need to extract all the string before the dot, excluding .CEL
Just an example: for the first row, I need to extract `11:00223`, for the second row `9399` and so on...

I tried using the library `stringr` and then the following command line:

str_extract(df$IID, "[[:digit:]]+")

but `[[:digit:]]+` extracts only number before the dot, while I need all the string before the dot, included letters of the alphabet.

Any idea about how to extract all the string before the dot?

Thank you!

Answer: gsub("\\..*", "", string)
Question: Hi,

I just started to work with single end reads, which are already trimmed for adapter sequences and quality.
Do I have to trimm the reads now to the same length of e.g. 100nt for mapping them with STAR?
Is there a negative effect, if I don't?

Answer: If the qualities are ok and there are no adapters you can proceed with mapping.
There is a recent paper about trimming of RNAseq data and its possible consequence on downstream analysis - https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4766705/
Question: I just do some practice on ballgown and stringtie, and I got some GTF file and ballgown`s file.

However, I just find I can't use R or something else to deliver a excel file which contain FPKM.

I think the excel file I want maybe looks like this:

gene_id FPKM
A 124
B 541
C 122

Please help me

Thanks a lot :)

Answer: You can use standard unix commands for this as follows:

$ cat stringtie.txt
1 StringTie transcript 337772 338047 1000 + . gene_id "Zm00001d027250"; transcript_id "Zm00001d027250_T001"; cov "66.076088"; FPKM "8.407302"; TPM "14.141658";
1 StringTie exon 337772 338047 1000 + . gene_id "Zm00001d027250"; transcript_id "Zm00001d027250_T001"; exon_number "1"; cov "66.076088";
1 StringTie transcript 426764 432130 1000 + . gene_id "Zm00001d027254"; transcript_id "Zm00001d027254_T001"; cov "2.043083"; FPKM "0.259955"; TPM "0.437262";
1 StringTie exon 426764 426798 1000 + . gene_id "Zm00001d027254"; transcript_id "Zm00001d027254_T001"; exon_number "1"; cov "0.000000";
1 StringTie exon 426869 426970 1000 + . gene_id "Zm00001d027254"; transcript_id "Zm00001d027254_T001"
$ grep 'FPKM' temp.txt | cut -f9 | sed -r 's/gene_id "([^"]*).*FPKM "([^"]*).*/\1\t\2/g' | sed '1i#gene_id\tFPKM'
#gene_id FPKM
Zm00001d027250 8.407302
Zm00001d027254 0.259955
Question: Hi all,

I search the coordinates of PAR3 on chrX and on chrY. Wikipedia (please see link attached) tells me, that the genes involved are PCDH11X and TGIF2LX on X (@ ~ 90MB); PCDH11Y and TGIF2LY on Y (@ ~ 5 MB). The coordinates of these gene serve me as a rough orientation, but I search an 'official definition' for example from the UCSC or any other official instituion.

Do you know of any?

Thx,

Best,
Mathias.

Wikipedia: http://en.wikipedia.org/wiki/Pseudoautosomal_region

Answer: I did not know PAR3 until now. The claim of PAR3 is made by one paper last year which has only been cited once. In the paper, the author is calling XTR (X transposed region) as PAR3. I have not read through the paper, but I think the evidence described in the abstract is very weak. Correlated copy number changes are not sufficient to define PARs - ordinary segmental duplications are also associated with non-homologous recombinations sometimes. Biologically, a PAR in the middle of chrX is also weird. If there is a single crossover between chrX and chrY in PAR3, would that create a huge hybrid chrX-chrY chromosome? I would take PAR3 with a pile of salt.

On the other hand, XTR is a well-known segmental duplication from chrX to chrY. You can find its approximate boundaries from UCSC self-chain:

https://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&position=chrX%3A87015838-93805574&hgsid=386334885_AcoZWTXT6wlO6yJPaqfacRL0k0Cl

Note that there are massive inversions and deletions between the chrX/chrY copies (there are papers about it, but I cannot find them now). It is not possible to precisely define the exact boundaries.

EDIT: XTR evolution: <http://www.nature.com/nature/journal/v423/n6942/extref/01722/figures/fig5.pdf>
Question: Hello everyone,

My question comes in 2 part.

First, I have a R data frame of 1.3M human SNPs (all with a RS number) and would like to:

- Update their genomic position (the position I have is based on an old
hapmap reference, it would be nice to have it with Hg38 coordinates)
- Annotate their effect (as in the CONTEXT column in the GWAS catalog, including for example : 3_prime_UTR_variant 5_prime_UTR_variant, downstream_gene_variant, intergenic_variant, intron_variant, missense_variant, non_coding_transcript_exon_variant , regulatory_region_variant and a few other categories)

I looked at many R packages, including rsnps and BSgenome, but none of them was able to extract dbSNP information (only rsnps could extract their position but on a very small subset). I am aware SnpEff could (maybe, not sure) work with rsID but I would like to stick with R for that part.

----

Second question, I am using SnpEff on a vcf-type file of SNP list, and I would like to get the annotation on only the coding variants and for them whether if it is Synonymous or Missense/LoF for the non synonymous.

I ran SnpEff and got a shitload of extra information and it's hard to make sense of it. So I used the filters -no-downstream -no-intergenic -no-intron -no-upstream -no-utr -no EffectType (low).
Fore some reason, the output still contained a lot of the stuff I wasn't interested about, apparently the filters weren't the good ones.

I anybody can help me with those questions, that would be great! I am starting to run a bit dry on the answers Google can tell me.

Many thanks,

Answer: R code:

rsid=c("rs123","rs150")
library(biomaRt)
ensembl_snp=useMart("ENSEMBL_MART_SNP", dataset="hsapiens_snp")
getBM(attributes=c("refsnp_source",'refsnp_id','chr_name','chrom_start','chrom_end',"consequence_type_tv", "clinical_significance"), filters = 'snp_filter', values = rsid,mart = ensembl_snp)

output:

> getBM(attributes=c("refsnp_source",'refsnp_id','chr_name','chrom_start','chrom_end',"consequence_type_tv", "clinical_significance"), filters = 'snp_filter', values = rsid,mart = ensembl_snp)

refsnp_source refsnp_id chr_name chrom_start chrom_end consequence_type_tv
1 dbSNP rs123 7 24926827 24926827 intron_variant
2 dbSNP rs150 7 24971033 24971033 intron_variant
clinical_significance
1 NA
2 NA

coordinates seem to match with dbSNP records.
Question: Hi, I'm an undergraduate student. Please help me.

I want mapping-bam-file from sra-dataset from NCBI in order to analyze heterogeneity of mouse ESC.

The dataset is generated from a paper below.(GSE60749)

*Roshan M.Kumar et al. Deconstructing transcriptional heterogeneity in pluripotent stem cells. Nature(2014)*

Please teach me how to get adapter sequence used in this experimentation

and what I should use for quality control(Prinseq?ShortRead?).

For example, I want to try GEO Sample [GSM1486817][1] sra file.

Can anyone give me process of quality control of this sra?

I thank you for reading it through.

Any help will be appreciated.

[1]: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM1486817

Answer: 1. Download data
2. Convert from` .sra `format to .`fastq` with **SRA Toolkit**: http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software
3. Check the quality using **fastQC** package: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/. In the generated report, you will see if your sequences have adapters and if so, their names.
4. Remove the adapters and bad quality sequences using:
- **Trimmomatic**: http://www.usadellab.org/cms/?page=trimmomatic
- **AdapterRemoval**: https://github.com/mikkelschubert/adapterremoval
- **Cutadapt**: http://journal.embnet.org/index.php/embnetjournal/article/view/200/479
- ... and many other tools. I would suggest you Trimmomatic.

These are the general steps that you should follow in order to perform a quality control of the raw data.

Hope it helps.
Question: Hi. I'm trying to do random sampling with R.
When I run

library(GenomicAlignments)
library(rtracklayer)

## your bam file
bam.file <- 'your_bam_file.bam'

## is it paired end or not
paired <- TRUE

## read in the bam file
if(paired){
bam.input <- readGAlignmentPairs(bam.file)
}else{
bam.input <- readGAlignments(bam.file)
}

## select your number of subsets
num.samples <- 100

## subsample WITH replacement
bam.subsample <- sample(bam.input,num.samples,replace=TRUE)

## export new bam file
export(bam.subsample,BamFile('bam_subsampled.bam'))

I get the following error message

[E::hts_idx_push] Region 554468587..554468738 cannot be stored in a bai index. Try using a csi index with min_shift = 14, n_lvls >= 6
value[[3L]](cond) でエラー: 'asBam' failed to build index
file: bam_subsampled.bam
SAM file: 'bam_subsampled.sam'
追加情報: 警告メッセージ:
.make_GAlignmentPairs_from_GAlignments(gal, strandMode = strandMode, で:
79328 alignments with ambiguous pairing were dumped.
Use 'getDumpedAlignments()' to retrieve them from the dump environment.

I searched for csi option in Rsamtools or GenomicAlignmets, but I wasn't able to find that.

How can I solve this error ??
Thanks.

Answer: Simply add `index = FALSE` to avoid building the index.

`export(bam.subsample,BamFile('bam_subsampled.bam'), index = FALSE)`

If you need the index you can create it with samtools at any time (`samtools index -c bam_subsampled.bam`)

Also asked for some details over at the BioC forum, see what they say. https://support.bioconductor.org/p/126603/
Question: Hi, I have been able to create my own linear workflows (A -> B -> C -> D) with Snakemake. However, now I would like to include some optional steps (X) that should be only executed if the user specifies it. Briefly, most of the times C will take as input the output from B, but sometimes I would need X to take as input the output from B, and then C to take the output from X.

![enter image description here][1]

After looking at the documentation I have not been able to figure out how to do this. I don't even know if this is feasible, or if there is other approach that fits better. I would appreciate some guidance here.

Thanks!

[1]: https://i.ibb.co/3TkqHHx/Diagram.png

Answer: The input to `C` can be a function that defines files based on the wildcard for the output from `C`.
You could write a function that decides what the input to C should therefore be.

When you say that "sometimes I would need X to take as input the output from B and then C to take the output from X", do you mean that the optional use of X is decided for a given sample within your workflow (sample 1 might pass through X, but sample 2 might not need to), or that the whole workflow should optionally use X based on some config/argument (for a given experiment, choose to run all samples through X)

Something like this

# optionally use X for every sample passing through the workflow
def input_for_c(wildcards):
# requires a config containing switches for the whole workflow
if config["Use X"]:
return "./data/X/{}".format(wildcards["sample_id"])
else:
return "./data/B/{}".format(wildcards["sample_id"])

# optionally use X for the current sample
def input_for_c(wildcards):
# requires a `sample_config` containing switches for each separate sample
sample_id = wildcards["sample_id"]
if sample_config[sample_id]["Use X"]:
return "./data/X/{}".format(wildcards.sample_id)
else:
return "./data/B/{}".format(wildcards.sample_id)

rule all:
input: expand("./data/C/{sample_id}", sample_id = SAMPLES)

rule B:
output: "./data/B/{sample_id}"
...

rule X:
output: "./data/X/{sample_id}"
...

rule C:
input: input_for_c
output: "./data/C/{sample_id}"
...
Question: What is the difference between Phred and GQ(Genotype quality)?
Which are the values for a good cut-off?

Answer: Phred scale is usually used for representing the base error probability in lograrithm format:

$Phred score = -10 log_10 error probability$

genotype quality is how certian the genotype calling algorithm is about the genotype called. It is usually displayed in 10*log(relative_score_of_genotype) format too.
Question: Hi,everybody,

I find that the lastest version of gene in NCBI is GRCh38,I could find GRCh37 for on-line browser version. But I can not find the download version.In the download page, The only version is GRCh38. Anyone know where to download GRCh37 download files in NCBI?

Answer: under ftp://ftp.ncbi.nlm.nih.gov/genomes/Homo_sapiens/ARCHIVE/BUILD.37.3
Question: Dear Biostar Community,

What is the best tool for read generation from given reference sequences or any that you use/used and it works/worked as expected?

The properties of desired tool:

- reads generate of given insert size and read length with uniform coverage
- handle uniform coverage on circle genomes
- handle large genomes without crashing
- generation of Illumina reads
- generation of PacBio reads
- generation of reads without errors
- generation of reads contaminated by technical sequences
- generation of read-through adapters artefacts

I used ART (http://www.niehs.nih.gov/research/resources/software/biostatistics/art/) but it looks like it uses too old model for Illumina.

Answer: [randomreads][1] from BBMap.
PacBio reads can be generated [this][2] way.

[1]: http://seqanswers.com/forums/showthread.php?t=48988
[2]: https://sourceforge.net/p/bbmap/discussion/general/thread/25863b5d/#738f/91e5
Question: Hi,

I have one question, I need to get reads after filtration. Than normally I get health reads (health = reads after filtration, satisfy the conditions), but I need to get reads, which I remove.

For example I set filtration to get reads with mapQ bigger than 10, I used `samtools view -bg 10 ...` and I get reads with mapQ bigger than 10, but I want reads with mapq lower than ten.

Could you please help if it is possible. And other option, compare the raw bam file and filtrate bam file and get the difference between them.

Could you help me pls or do you know some tool to do this?

Answer: Just use awk or another scripting language, for example `samtools view foo.bam | awk '{if($5<10) print $0}'`. You can make that a bit more complicated if you want to keep the header (useful for piping back into samtools for conversion to a BAM file), but that's enough to give you the idea.

BTW, pysam is useful in general for more complicated tasks.
Question: I have a vcf file (format VCFv4.0), generated by GATK pipeline starting from Illumina reads.

I need to convert it to 23andme file format. Example of the 23andme format:

# rsid chromosome position genotype

rs4477212 1 82154 TT

rs3094315 1 752566 TC

rs3131972 1 752721 AA

rs12124819 1 776546 AC

I am having problems with plink2 `--recode 23 cannot be used with multi-char alleles`. Plink was recommended earlier here https://www.biostars.org/p/198034/#274523

I tried then to modify the vcf to remove multi-char alleles using [VcfMultiToOneAllele,][1] which did a great job but the output file, even though it looks like a vcf, it was not recognised as such by plink2 `no genotype data in .vcf file`. Any other tool up to the task?

Thanks for any help.

[1]: https://lindenb.github.io/jvarkit/VcfMultiToOneAllele.html

Answer: OK, it seems I have solved it using:

plink2 --vcf [vcf file] --snps-only --recode 23
now thinking how to include single point deletions and insertions to the output file, because those are missing
Question: I have a data frame of omics data. Gene ids in rows (931), and samples in columns (15).

> dim(my_data) # (rows columns)
[1] 931 16

I created heatmap using `library(gplots)`

cn=colnames(gdf1)[c(13:15,1:12)]

col <- colorRampPalette(c("red","yellow","darkgreen"))(30)

heatmap.2(as.matrix(gdf1[,cn]),
dendrogram = "row",
Colv = FALSE,
Rowv = TRUE,
scale = "none",
col = col,
key = TRUE,
density.info = "none",
key.title = NA,
key.xlab = "Abundance",
trace = "none",
margins = c(7, 15))

However, in the heatmap, I see only few genes. Since it has 900 ish genes. How can I export what are the clustered genes in each cluster in the same order as in the heatmap?

Also, how can I reduce the size of colorkey?

Thank you.

Answer: **Edit 13th September, 2019:**
To additionally see how to extract clusters of genes from the heatmap dendrogram, zoom down to this later comment: https://www.biostars.org/p/398548/#398618

---------------
---------------

Hello,

First, create random data
----

mat <- matrix(rexp(200, rate=.1), ncol=20)
rownames(mat) <- paste0('gene',1:nrow(mat))
colnames(mat) <- paste0('sample',1:ncol(mat))
mat[1:5,1:5]
sample1 sample2 sample3 sample4 sample5
gene1 0.6247039 3.020142 8.303563 6.482744 0.59547154
gene2 2.6650871 3.375123 5.778222 19.410709 0.07966728
gene3 4.6343755 5.491166 8.716883 9.490372 29.03157875
gene4 13.6086878 3.632815 10.688699 1.263853 2.54216953
gene5 2.4060078 14.283380 8.592085 3.998141 0.25853135

Generate a heatmap and save it to `out`
---

out <- heatmap.2(mat)

Obtain list of genes, ordered as per heatmap (from bottom, up):
---

rownames(mat)[out$rowInd]
[1] "gene2" "gene9" "gene7" "gene8" "gene4" "gene5" "gene3" "gene6"
[9] "gene1" "gene10"

Plot the row dendrogram on its own:
------

plot(out$rowDendrogram)

Change colour key size
-------

Use `keysize` parameter

----------------------------------
-----------------

See also here for pheatmap: https://www.biostars.org/p/287512/#287518

Kevin
Question: Hello,

I have a set of 5kb sequences upstream of a gene from different primates, and I would like to know what motifs are in enhancers and promoters of primates vs the 5kb upstream sequences of a set of the same gene but for rodents.

Would the analysis with MEME suite be correct with the "Discriminative mode", and putting the rodent sequences in control sequences and the primate sequences in the primary sequences? To contrast, since the normal method looks for similarities.

If not, how would this analysis be carried out?

Thanks!

Answer: I think the discriminative mode is more to look for enrichment of the same motif. But different species may have the same motif with different information values.

So one thought is, given sets of orthologous genes, you might perhaps perform a MEME scan for conserved motifs over their respective upstream promoter sequences. Pairs of per-species MEME PWMs could be compared with TOMTOM to query for significantly similar motifs. Putting results into a grid could help with comparing all species pairings at a glance.

Another thought is to use existing (MEME-formatted) PWMs with FIMO to generate whole-genome scans for published motifs, comparing hits over ortholog promoter regions.
Question: Hi,

I have a set of samples and log10 total read counts of them in a data frame like below

> head(a)
log10_total_counts
A1 6.468503
A10 6.565213
A11 6.752139
A12 5.078598
A2 6.277342
A3 6.473411
>

![enter image description here][1]

[1]: https://i.ibb.co/P6mmKc5/Picture1.png

For instance red dots have less than 1000000 counts and black dots have more thsn 1000000 counts.

How could I repreduce this plot please?

Answer: another plot in ggplot2 with simulated data:

df=data.frame(genes=paste0("gene_", seq_along(1:75)), reads=round(rnorm(75, 110,2),0))
library(ggplot2)

ggplot(df,aes(genes, reads, color = ifelse( reads < 110, "Fail", "Pass")))+
geom_point() +
scale_y_continuous(limits = c(0, 200))+
scale_color_manual(name="QC", values = c("red","darkgreen"))+
geom_hline(yintercept = 110, color="red")+
theme_bw()+
xlab("Sample")+
labs(y=expression("log"[10]("Total reads")))+
theme(axis.text.x = element_text(angle = 45,hjust = 1),
legend.position = "bottom")
Question: I have a list of >100 gene names and want to get the `locus_tag`s. Is there a way I can do this using its genbank file and Biopython?

For example, the genbank file has

```
/gene="murE"
/locus_tag="BSUW23_07815"
```

If my list just has `murE`, I'd like it to print out the corresponding `BSUW23_07815`

Answer: Something like this should work:

```
from Bio import SeqIO
​genbank_file = "example.gbk" # insert your filename here
wanted = ["murE", ...] # or load all your 100 genes from a file
for record in SeqIO.parse(genbank_file, "genbank"):
for f in record.features:
if f.type == "CDS" and "gene" in f.qualifiers:
gene = f.qualifiers["gene"][0]
if gene in wanted:
print f.qualifiers["gene"][0], f.qualifiers["locus_tag"][0]
```

See also http://www.warwick.ac.uk/go/peter_cock/python/genbank/
Question: Hello all,

I have about ~150 sub directories in a directory. Each of those subdirectories contains multiple fasta files. I want to rename all the fasta headers with adding the name of the file.

I tried using this command:

for i in ./*mg/*.faa; do gawk -i inplace '/>/{gsub(">","&"FILENAME"_");gsub(/\.faa/,x)}1' $i; done

This command however renames the fasta headers with sub-folder name (*mg) and the filename. How do I modify this command to just include the filename??

Answer: find /*mg -type f -name "*.faa" | while read F; do sed "/^>/s|\$| $(basename $F)|" "${F}" > "${F}.new" ; done
Question: Hello, I have a text file likes this:

GTCAAAGATCAGATGGTTGTAGATGTGTGGTGTTATTTCTGAAGCCTCTGTTCAAGGAAGGCAGTTTCTTATGAA

I want to convert it to a fasta file. Can we just manually add a ">" in the fist line?

Thanks

Answer: A good explaination of the FASTA format is here:

http://en.wikipedia.org/wiki/FASTA_format

Some applications will use the word immediately after the ">" symbols as an "ID" so make sure it means something to you. You that you can put multiple sequences in the one file, just start each new one with a ">" line, and use a unique ID for each one.
Question: In the past when I've ran Abyss, I've given it multiple sequencing libraries from multiple input fastq.gz files (~36 separate files).

In this current run, I've concatenated all my sequencing libraries together into two huge fastq.gz files (one for each of the pair). The main reason was that I performed digital normalization and the output was two huge concatenated files.

The issue now is that I am finding Abyss to be taking a lot longer to finish assembling these normalized reads (about 1/3 of the original size).

From reading the abyss logs, I noticed in the past there was parallel reading of multiple input fastq.gz files. Is it possible that since I just have two huge files now, I am not taking advantage of this parallel reading? And a significant amount of time is being used just to read in files?

In my past runs with multiple inputs, Abyss would finish reading in the reads in 6-8 hours (I determine this by check when the pruning started to get logged). It has been almost a full day now with the current run and it still haven't finished reading in the reads.

I am running Abyss on an ad hoc AWS cluster (starcluster) with 20 nodes. In past runs with multiple input files, I've used abyss 1.52. In this current run with two huge files, I am using abyss 1.9. Could it be a difference between versions also?

Answer: > From reading the abyss logs, I noticed in the past there was parallel reading of multiple input fastq.gz files. Is it possible that since I just have two huge files now, I am not taking advantage of this parallel reading? And a significant amount of time is being used just to read in files?

Yes. This issue is on our radar. Feel free to open a feature request issue on GitHub.

> In past runs with multiple input files, I've used abyss 1.52. In this current run with two huge files, I am using abyss 1.9. Could it be a difference between versions also?

No.

Cheers,
Shaun