cannin/biostars_qa huggingface dataset:
Question: Hi guys, I have thousands of fasta files each with one sequence only. I need to run a command in my linux/unix system to perform certain job and loop that command through each file and get the output in a single file. Is there a way to do that? I was suggested to use bash loop. For example, if I need to translate a DNA sequence in a fasta file and if the program (sixpack in this case) translates only one sequence at a time, how do I do this for multiple submission of files? I want to avoid doing this individually:

sixpack <fasta1.fasta
sixpack <fasta2.fasta
sixpack <fasta3.fasta

Instead, I want the bash(or perl/python script) to follow the same command looping through each fasta file. Thank you for your help!

Answer: for a larger/longer process that needs to be parallelized,
use GNU parallel: http://www.gnu.org/software/parallel/

parallel sixpack '<' {} ::: *.fasta

or GNU Make with option -j
Question: I have data frame like this about 250 gene i taking a small subset of it

gene `WT` `Treated`
ACTL6A	4.5167211886	4.4304190975
ACTL6B	-4.3730103378 -9.800779238
ACTR5	1.3871253811	1.2086608355

I am plotting heatmap since I m using pheatmap since plotting requires a matrix that is I have to use the only the numeric matrix
I m doing something like this
`rownames(VD3) <- data[,1]`

`mat <- as.matrix(data[,-1])`

I can plot the heatmap but Im missing my gene name labels , How can i include it in my final plot ,any suggestion ?And Im using the pheatmap library

Answer: Your code is a bit confusing.
You can assign the rownames to the matrix.

Mat = as.matrix( data[,-1])
rownames(Mat) = as.character( data[,1])
head(Mat)
class(Mat)

Did this help you?
Question: Hi there,

this is my first post at BioStars. I am a new starter bioinformatician. I have a problem that should be easy to solve but I cannot sort out, so I would be glad if you could help me with this issue.

The question is easy. I have a VCF file with genotype data of many samples. It contains SNPs in the rows and the columns are of a typical VCF (#CHROM,..., INFO) followed by the Ids of the samples. I would like to filter out the non-European samples according to these genotype data, using HapMap. I was told I had to do a PCA. I have tried several tools for this, Shellfish, Beagle, SNPRelate, but I could not solve the problem. With SNPRelate, I could do the PCA, but this just clusters samples that are unlabeled and I need to associate them to HapMap populations (CEU, YRI, JPT, CHB). On the other hand, Shellfish returns me a non-informative error when it is running:

Exception: command gtool -P --ped shellfish-temp-15479/146134504516.ped --map shellfish-temp-15479/146134504516.map --og shellfish-temp-15479/146134504516.gen --os shellfish-temp-15479/146134504516.sample --discrete_phenotype 0 >> shellfish.log exited with code 256 (1)

And in file shellfish.log:

...
Note: No phenotypes present.
--recode to plink.ped + plink.map ... done.
Unknown parameter: 0

What steps and tools would you recommend to follow? I can use any tool you think is suitable for this.

Sorry for this, it may be an easy problem, but I have spent 2 days trying several tools.

Many thanks.

Answer: Hi,

Dont have experience with those tools but I usually use plink.

First convert your vcf file to plink file using for example vcftools (http://vcftools.sourceforge.net/)

./vcftools --vcf input_data.vcf --plink --out output_in_plink

Then download the hapmap data from plink: http://pngu.mgh.harvard.edu/~purcell/plink/dist/hapmap_r23a.zip

Then extract snps lists from both datasets and filter based on the these snplist to get only overlapping snps.

plink --bfile fileA--write-snplist --out list1 --noweb
plink --bfile fileB --extract list1 --noweb --out fileB_filtered --make-bed

Then merge the files and make a mds plot

plink --bfile fileA_filtered --bmerge fileB_filtered.bed fileB_filtered.bim fileB_filtered.fam --noweb --out merged --make-bed

plink --bfile merged --mds-plot 2 --noweb --out mds

This can be easily plotted using R or even Excel... and then you can see which samples are derived from which ethic background.
Question: In the new version of FastQC (0.11.2) I cannot find the `fastqc_data.txt` file. Has anyone an idea of how to create it? Is there an option to set? Or is it gone for good?

Answer: Problem solved. I added `--extract` to automatically extract the .zip file and then the `fastqc_data.txt` file appears again.
Question: Dear all,

I am working on vcf data with *bcftools*. It is a big data set and I would like to see the list of samples are including in this vcf file, what's the easy way with bcftools or vcftools?

Thanks and with best regards

Xinhui

Answer: I got it:

bcftools query -l input.vcf

From [manuals - query][1]:
> `-l`, `--list-samples`: list sample names and exit

[1]: http://www.htslib.org/doc/bcftools.html#query
Question: Is there any way to batch download UniProt proteomes based on a list of proteome IDs (e.g. UP000011712).

I have a list of about 1,000 proteome IDs from UniProt. I know it is possible to do this with taxon IDs but I do not have these.

So, is it possible to batch download proteomes based on proteome IDs, or is there any way to convert proteome ID to taxon identifier?

Thanks for any help.

Answer: You can do this programmatically using URLs like
https://www.uniprot.org/uniprot/?query=proteome:UP000005640&format=fasta (FASTA format)

Use format=txt, format=xml etc for other download formats.

There is a related FAQ for programmatic access with a code example to download the UniProt reference proteomes for all organisms below a given taxonomy node in compressed FASTA format: https://www.uniprot.org/help/api_downloading
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: GNU parallel and seqkit answer for posterity.

```
parallel 'seqkit replace -p "(.+)" -r "\$1 {/}" {} > {= s/\.faa/.renamed.faa/ =}' ::: $(find ./*mg -name "*.faa")
```
Question: I am trying to identify genomic windows in which the sequence is enriched in CpHpG or CpHpH motifs. Where H is any base that is not a G [ATC].

I have done this already for the di-nucleotide motif CpG by setting a threshold value for Obs/Exp as per EMBOSS's cpgplot, wherein:

Observed = Count of CpG di-nucleotides in a window
Expected = (number of C's * number of G's) / window length​

It is easy enough to count the number of CpHpG motifs in a given window, as these cannot overlap. However, I'm not sure how to best calculate the number of CpHpG's that would be observed by chance in a sequence of known length where we know the proportion of all bases.

Any suggestions appreciated!

Answer: Hi Adam,

A while ago I wrote a java program to test whether a regular expression appears in a sequence more often than expected by chance given the nucleotide composition of the sequence. P-value for enrichment is obtained by shuffling the sequence either by simple random sampling with replacement or by using a Markov model of given order.

I put it [here](http://bioinformatics-misc.googlecode.com/svn/trunk/eclipse/RegexEnrichment/RegexEnrichment.jar), source is [here](http://bioinformatics-misc.googlecode.com/svn/trunk/eclipse/RegexEnrichment/src/regexEnrichment/).

It might suite the task is addressing.

Show the help:

java -jar RegexEnrichment.jar -h

A simple test fasta sequence:

echo ">seq1" > test.fa
echo "AAAAAAAAAAAACTGAAAAAAAAACAGAAAAAAAAAAAAACAGAAAAAAAAAAAACTGAAAAAAAAACAGAAAAAAAAAAAAACAG" >> test.fa

Example:

Get the enrichment of the CHG pattern in windows of 40bp. Shuffle sequences using a first order Markov model.
By default, it returns a line in bed format for each pattern found. If you are interested only in windows, pipe the output to awk & uniq:

java -jar RegexEnrichment.jar -w 40 --norevcomp -omc 1 -p 'C[ACT]G' -f test.fa \
| awk 'BEGIN{OFS="\t"}{print $1, $7, $8, $5}'\
| uniq

seq1 0 40 0.3750
seq1 40 80 0.0800
seq1 80 86 0.2200

Last column is an empirical p-value for the pattern being enriched by chance in each tested window.
(Windows without pattern are not outputted).

Any feedback welcome...
Question: This is inspired by a comment on [this post](https://www.biostars.org/p/9508530/#9508566) by /u/dariober:

>Basically, you are asking snakemake to produce one log file per SRR id and this log file is produced by rule download_files. As a side product, download_files will give you the actual fastq files (things could be done differently in this respect but hopefully this will help...)

I have a list of SRA files that I need to download to split reads. The goal is to get from `[SRR1, .....]` a list of files `SRR1_1.fastq`, `SRR1_2.fastq` in a predetermined folder. To that extent, I wrote the following file:

SRA_MAPPING = read_dictionary()
SRAFILES = list(SRA_MAPPING.keys())[1:]

RawSampleFolderName="raw_samples"
RawSampleFolder=RawSampleFolderName+"/"

rule all:
input:
expand("{RawSampleFolder}{srafiles}_1.fastq",srafiles=SRAFILES, RawSampleFolder=RawSampleFolder),
expand("{RawSampleFolder}{srafiles}_2.fastq",srafiles=SRAFILES, RawSampleFolder=RawSampleFolder)

rule download_srafiles:
output:
expand("{RawSampleFolder}{srafiles}_1.fastq",srafiles=SRAFILES, RawSampleFolder=RawSampleFolder),
expand("{RawSampleFolder}{srafiles}_2.fastq",srafiles=SRAFILES, RawSampleFolder=RawSampleFolder)
params:
download_folder = RawSampleFolderName
shell:
"fasterq-dump {wildcards.srafiles} -O {params.download_folder}"

(This is a proof of concept, I'll dump the global variables in `config` as soon as I can). The nutshell is that I have a list of SRA files read from a list, and I have a preset download folder. I use `fasterq-dump` to get split read files - from `SR1` to `SR_1.fastq` and `SR_2.fastq`, In a snakemake fashion, I'd like for the rule `download_srafiles` to have `output` to be fastq files. My previous solution was given in the linked posted, but said solution made a `.log` file as an output and retrieved files as a side effect - I'd like to skip the part where I need `log` file. Since all I have is a python list, I skip the input. Running the above file does not download samples. Instead I get the error:

'Wildcards' object has no attribute 'srafiles'

So what is it that I'm doing wrong?

Answer: The reason you are getting this error message is that you have no wildcards in your `download_srafiles` rule. When you use `expand()`, the function simply replaces the values enclosed in {}'s, so no wildcards remain. This is a bit confusing, but keep in mind that {}'s within expand means "variable to replace", whereas in normal inputs/outputs {}'s indicate wildcards.
In fact, you shouldn't use expand in your second rule - this is already taken care of by the rule `all`. So the code for the second rule should look something like:

rule download_srafiles:
output:
"%s{srafiles}_1.fastq" % RawSampleFolder,
"%s{srafiles}_2.fastq" % RawSampleFolder
params:
download_folder = RawSampleFolderName
shell:
"fasterq-dump {wildcards.srafiles} -O {params.download_folder}"

Finally, I am not familiar with fasterq-dump, but I highly recommend [Kingfisher][1], which can also download from ENA which is much faster.

[1]: https://github.com/wwood/kingfisher-download
Question: Hi,

I want to ask in analysis of rnaseq data to detect differential expressed genes, etc., can technical replicates and biological replicates put together to compare?

Thanks,

Ch

Answer: The standard recommendation is to combine your technical replicates into one "lane" by simply adding their reads together. Use your biological replicates to estimate variance/dispersion.
Question: How to derive pan proteome of 28 microbial species?

I already have core genes / orthologs between all the species.

I also have unique genes of all species.

I only had problem with getting accessory genes/ proteins.

I am not willing to do blast 784 times.

help pls

Answer: First, pan-proteome is not a widely-used term, nor is it very well defined. I assume you mean it in the sense used in this article:

core (found in all) + accessory (found in > 2) + unique (found in 1)

I don't think the 3 sets together are interesting or useful ; pan-proteome is just a term used to describe those 3 sets.

Second, why are you "not willing" to run 784 BLAST searches? Perhaps you mean "not able"? If you're able but not willing, perhaps consider another career ;)

Third, BLAST is not necessarily the best tool. You may find it more useful to cluster the protein sequences using e.g. CD-HIT. In fact, the output from that may come close to giving you the groupings that you need.
Question: Hi every one i am doing some plamid genome assembly with spades. After assembly i used SSPACE for scaffolding. But there are some gap each of the draft genome. I can fill the gap by PCR. But i want to reduce the number of gap insilico?? Can anyone suggest how to reduce the gap ?? If i do mapping the reads with contigs will it be give any promising result???

Answer: If you use paired reads you can take a look to: GapFiller

I haven't test it myself, though.
Question: Hello all,

Using time series **(two condition with 4 time point)** data in DESeq2, I got warning while fitting by GLM ??

bckCDS <- DESeqDataSetFromMatrix(countData = x, colData=ExpDesign, design= ~condition.label + condition.Time + condition.label:condition.Time)

ddsLRT= DESeq (bckCDS, test =LRT, reduced = ~condition.label + condition.Time)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
NOTE: fitType=parametric, but the dispersion trend was not well captured by the
function: y = a/x + b, and a local regression fit was automatically substituted.
specify fitType=local or mean to avoid this message next time.
final dispersion estimates
fitting model and testing

May I have fit it by "local", indeed its work by this however, my question is why its not working by `fitType='parametric'`. What is advantage and disadvantage of it?

Thanks

Answer: Depending on the version you are using, this used to be a warning but is now just a message/note.

The printed message here is literally all you need to know. The software cannot find parameters a and b to fit the dispersion over mean trend using that formula (read the manuscript for more information about the point of modeling dispersion over mean).

So a local regression was automatically substituted. There is no further action required by the user.
Question: I have been working on a script for identifying frequent 'end points' of paired end reads from a SAM file on a genome for determining uneven cut frequencies, which I need to integrate the orientation and 1st/2nd read information into. The 2nd column of a SAM file gives the detail on the orientation with a number (as described [here][1] and [here][2]).

I hadn't used a bioperl import method because all I was doing was pulling the genomic start position of the read and building a hash by frequency, but I need to now tune it for directionality and mate pair.

My question is whether there is a better way to do this than:

```
## Numbers are random examples, not correctly done at the moment.
if ($line[1] == 99){
$coordinate = $start
}elsif($line[1] == 147){
$coordinate = $start - $readlength
}elsif($line[1] == 83){
$coordinate = $start + $readlength
}.........etc
```

Given that a picture paints a thousand words, basically what I want is to pull out the coordinates of the red circles from my SAM file. It's part of a bigger perl script, so I'd prefer to avoid other tools, but if they exist I should look at them.

[][3]

[1]: http://bowtie-bio.sourceforge.net/manual.shtml#sam-bowtie-output
[2]: http://http://picard.sourceforge.net/explain-flags.html
[3]: http://i.imgur.com/EWT7mFA.png

Answer: You can parse the bitwise flags with the AND operator in perl. For example, the 5th bit in the bitwise flag denotes orientation:

```
if ($flag & 16) {
return 'negative strand';
} else {
return 'positive strand';
}
```

The 5th bit is 2^4 = 16. 1st bit is 2^0, 2nd bit is 2^1...
Question: Hi,

I'm analyzing **human cancer cell line** ATAC-seq data (PE 75bp Illumina) and have encountered a weird insert size distribution. For ATAC-seq PE data, insert size distribution is a good quality control of how well the experiment worked (This is how it would ideally look https://dbrg77.wordpress.com/2017/02/10/atac-seq-insert-size-plotting/). I have used several mapping strategies and always come to a similar distribution plot.

Sequencing is good quality with 50-60 million reads per sample with min. 3 replicates and mapping >80% for all samples. > 95% of forward and reverse reads are paired.

**Analysis**

java -jar $tool_path/trimmomatic-0.36.jar PE -phred33 -threads 20 $forward $reverse $forward_1_paired.fastq $forward_unpaired.fastq $reverse_paired.fastq $reverse_unpaired.fastq ILLUMINACLIP:$tool_path/adapters/TruSeq3-PE.fa:2:30:10 LEADING:10 TRAILING:10 MINLEN:50

bowtie2 --threads 20 --very-sensitive -x $bowtie2index -1 $forward_paired -2 $reverse_paired -S $filename_aligned.sam

**Plotting**

samtools view ATAC_f2q30_sorted.bam | awk '$9>0' | cut -f 9 | sort | uniq -c | sort -b -k2,2n | sed -e 's/^[ \t]*//' > fragment_length_count.txt

**Alternative way of plotting**

java -Xmx20g -jar picard.jar CollectInsertSizeMetrics I=$file O=insert_size_metrics.txt H=insert_size_histogram.pdf

[Insert size distribution][1]
*In the plot it should be written insert length, not fragment length.

Does anyone have an idea what went wrong here? Could this data be still useful to analysis? I was thinking to extract PE reads with insert size between 0-300nt and proceed only with these reads for downstream analysis?

[1]: https://ibb.co/frQb6U

Answer: For anyone looking at this post in the future, I will explain how I solved this issue:
The issue was in the raw data I was given. The data, as i mentioned above was PE sequencing and i was given two files per sample, marked as _1 and _2, and the information that these are forward and reverse reads. I processed them as such and gotten no errors and good mapping results. But in fact, what I realised later is that these were lane 1 and lane 2 files and both files contained both forward and reverse reads. Oddly i didn't get any errors while removing adapters or mapping, despite the fact that i treated them as if file 1 contained forward reads and file 2 reverse reads. The data now looks as it should, with the typical nucleosome phasing pattern and good quality.

So lesson learned, always double check the information you are given together with the files.
Question: I have a number of VCF files, where each VCF file possesses variant data for a single patient (this is the way Illumina provides their data). Is it possible to combine all of the data for the patients into one VCF file? If so, how? Can I use plink/seq to do this?!

Any suggestions and leads would be extremely helpful.

Answer: GATK *CombineVariants*, see:

- [CombineVariants - Combine variant records from different sources](https://www.broadinstitute.org/gatk/guide/tooldocs/org_broadinstitute_gatk_tools_walkers_variantutils_CombineVariants.php)

From the above link usage examples:

Merge two separate callsets

java -jar GenomeAnalysisTK.jar \
-T CombineVariants \
-R reference.fasta \
--variant input1.vcf \
--variant input2.vcf \
-o output.vcf \
-genotypeMergeOptions UNIQUIFY

Get the union of calls made on the same samples

java -jar GenomeAnalysisTK.jar \
-T CombineVariants \
-R reference.fasta \
--variant:foo input1.vcf \
--variant:bar input2.vcf \
-o output.vcf \
-genotypeMergeOptions PRIORITIZE \
-priority foo,bar
Question: I've got a list of codding genes and I would like to get the cellular location of the expressed proteins, all at once. I already looked up on uniprot but wasn't able to solve it.

Answer: http://compartments.jensenlab.org allows you to look up each gene individually, or download a flatfile with the location of each.
Question: Dear all,

Does STAR 2.6 `--quantMode GeneCounts` only quant unique mapping reads? How to specialize it only quant the unique mapped reads?

Thanks a lot.

Answer: Usa, STAR with `--quantMode GeneCounts` includes only uniquely mapped reads in the quantification.
Question: I wanted to know the best method/pipeline for gene annotation after draft genome assembly for gene annotation. Using BLAST is one way to go. Is there any tool (with good accuracy) that can annotate the genome??
Thank you in advance.

(My draft assembly is done using Masurca and SPades with paired-end reads, mate-pair reads and pacbio reads.)

Answer: I am assuming this is a Eukaryotic species. Depending on the specific species, you could use related-species' protein evidence as input along with a BUSCO-trained custom Augustus species and MAKER (http://www.yandell-lab.org/software/maker.html) to predict protein-coding genes. If there is a related, well-annotated species (ideally from ENSEMBL), you could use the Comparative Annotation Toolkit (CAT, https://github.com/ComparativeGenomicsToolkit/Comparative-Annotation-Toolkit). You could also run MAKER first and then align the **predicted cDNA transcripts** to your de novo genome and use that as input for Augustus PB in the CAT workflow.

edit: changed predicted proteins to cDNA transcripts
Question: I have a very basic question. I have created a local (mysql) database in my system. I would like to know what would be the best option for a front end. I had used html before and want to try out something new. I might also add some perl code in between for some computational work. Kindly help me in this regard..

Answer: Assuming that by "front end" you mean web interface, there are many possible solutions.

Since you mention Perl, you might want to look at a web framework such as:

- CGI::Application
- Catalyst

At the very least, you should be familiar with DBD::mysql; you can "roll your own" Perl CGI around that, perhaps using CGI.pm.