cannin/biostars_qa huggingface dataset:
Question: Hello,

There is a tool on the Galaxy platform which could potentially be very useful for the type of research that I am undertaking. The tool is 'Stitch MAF blocks' which can be found under the 'Fetch Alignments/Sequences' header.

The problem I have is that, because of the size of the data that I would be working with, I would need to try and run the function locally on my own system and using a command line (to enable access to local HPC resources)

I have installed the bx-python library (https://bitbucket.org/james_taylor/bx-python/wiki/Home) on my local system.

I think I am right in saying that these are the set of python scripts which essentially underlie the MAF manipulation tools (as well having many other functions) on the Galaxy platform.

However - the problem I have now - is that looking at all the files contained within the bx-python library - I can't easily discern a clear relationship between the python scripts contained within the library and the particular tool of interest on the Galaxy platform.

As a result, my question is the following:

Which Python files in the bx-python library correspond to the 'Stitch MAF blocks' function on the Galaxy platform?

Answer: The tool is executing the command defined here:

https://github.com/galaxyproject/galaxy/blob/master/tools/maf/interval_maf_to_merged_fasta.xml#L7-L9

so you are probably looking for file `interval_maf_to_merged_fasta.py`
Question: Dear all,

Could you help he please to modify my loop in bash to multiple alignment my FASTQ files?

I have a many fastq files R1 and R2 and I would like to align all those read in bash loop - lets say I have `X1_R1_001.fastq` + `X1_R2_001-fastq` and `Y1_R1_001.fastq` + `Y1_R2_001.fastq`. Sample X1 is pair-end read R1 + R2 and so on.

```
#!/bin/bash
for i in *fastq;
do tophat2 -o ${i%.fastq}tpout -G path/to/my/reference.gtf -p 8 path/to/my/bowtie_index ${i}R1_001.fastq ${i}R2_001.fastq --rg-id X1 --rg-sample X1 --rg-library rna-seq --rg-platform Illumina
done;
```

And I would like to also change my `--rg-id` tag with name of my fastq files (in this example `X1` in first loop and `X2` in second loop). Output folder should have the name of the sample too.

Please do you have any idea how to modify my bash loop?

Thank you so much for any ideas!

Paul

Answer: So, this is what I came up with, based on our discussion:

```
for i in $(ls *.fastq | rev | cut -c 13- | rev | uniq)
do
tophat2 -o ${i%.fastq}tpout -G path/to/my/reference.gtf -p 8 path/to/my/bowtie_index ${i}_R1_001.fastq ${i}_R2_001.fastq --rg-id ${i} --rg-sample ${i} --rg-library rna-seq --rg-platform Illumina
done
```

Hope this helps!

---

OK, this answer is gotten a little popular, and my shell knowledge has improved so I'm going to enhance this with some fail-safes to address commonly observed challenges.

First off, the `rev | cut -c 13-` is just asking for trouble. Let's substitute that with a `sed` so the command makes a little more sense. Here, I'm trying to remove all `_R[12]_001.fastq`.

`ls *.fastq | sed -r 's/_R[12]_001[.]fastq//' | uniq ` looks better. I've also edited the variable name from `i` so it makes a little more sense. And I've added some formatting changes plus quotes to avoid shell problems.

Here's the updated code:

```
for prefix in $(ls *R1*.fastq | sed -r 's/_R1_001[.]fastq//' | uniq)
do
tophat2 -o "${prefix%.fastq}tpout" \
-G path/to/my/reference.gtf -p 8 \
path/to/my/bowtie_index \
"${prefix}_R1_001.fastq" "${prefix}_R2_001.fastq" \
--rg-id "${prefix}" --rg-sample "${prefix}" --rg-library rna-seq --rg-platform Illumina
done
```

Note (21-Sep-2022): This will only work for one-pair-per-sample FASTQs. Multi-lane FASTQs will need to be concatenated so each sample has just one pair. The command can be modified easily for compressed FASTQ files.
Question: Dear Biostars,

I have a multisample vcf. I would like a to turn this into a per sample "alts" file as below.

The desired outcome would be :

chr10:12342:A:T SAMPLEID1:0/1 SAMPLEID2:1/1 ...
chr10:432234:G:C SAMPLEID4:0/1 SAMPLEID7:1/1 ...

I have tried :

bcftools query -f "%CHROM:%POS:%REF:%ALT\n[%SAMPLE:%GT\t]' i- 'GT="alt"' input.vcf

which gives me:

chr10:12342:A:T SAMPLEID1:0/1 SAMPLEID2:1/1 chr10:432234:G:C
SAMPLEID4:0/1 SAMPLEID7:1/1 ...

or

bcftools query -f "%CHROM:%POS:%REF:%ALT\t[%SAMPLE:%GT\n'] i- 'GT="alt"' input.vcf

chr10:12342:A:T SAMPLEID1:0/1
SAMPLEID2:1/1
chr10:432234:G:C SAMPLEID4:0/1
SAMPLEID7:1/1

or

bcftools query -f "%CHROM:%POS:%REF:%ALT\t[%SAMPLE:%GT\t'] i- 'GT="alt"' input.vcf

chr10:12342:A:T SAMPLEID1:0/1 SAMPLEID2:1/1 chr10:432234:G:C SAMPLEID4:0/1 SAMPLEID7:1/1

I can't quite get my head around the correct syntax and sure there is an easy fix! Your help would be much appreciated.

All the best

Answer: Figured it out!

bcftools query -f "%CHROM:%POS:%REF:%ALT[\t%SAMPLE:%GT]\n' i- 'GT="alt"' input.vcf

Just needed to move the \n out of the square brackets!
Question: I tried to download data from useast.ensemble.org, while the downloaded file only contains SNPs in chromosome 1. Is there a better way to do that?

Thank you in advance!

Answer: Ensembl's version at this link: https://ftp.ensembl.org/pub/release-92/variation/vcf/homo_sapiens/ They are split according to populations and also for 1000 genomes.
Question: I'm looking to do something like the following:

cwlVersion: v1.0
class: CommandLineTool

requirements:
- class: InlineJavascriptRequirement
- $import: Types.cwl

...

inputs:
an_input: string

outputs:
an_output:
type: 'Types.cwl#SpecialType'
outputSource:
valueFrom: ${
return {
'output_file': glob('the_tools_output.file'),
'this_should_be_associated_with_the_file_in_a_later_step': inputs.an_input
}
}

Is this possible? I see a possible solution here using an ExpressionTool

https://github.com/andersgs/cwl_flows/blob/master/limit_detection/seqtk_sample_PE.cwl

But I was wondering if it could be done inline in the tool output.

Answer: That should be possible, however the syntax for the output would look like this:

```
outputs:
an_output:
type: 'Types.cwl#SpecialType'
outputBinding:
glob: something
outputEval: |-
${
return {
'output_file': self,
'this_should_be_associated_with_the_file_in_a_later_step': inputs.an_input
}
}
```

In the expression for outputEval, `self` has the value of whatever is caught by `glob` pattern.
Question: I use picard software to mark duplicates, and here is my command:

```
java -d64 -server -XX:+UseParallelGC -XX:ParallelGCThreads=2 -Xms8g -Xmx16g -Djava.io.tmpdir=tmp -jar ./picard.jar MarkDuplicates I=input.bam O=out_markdup.bam METRICS_FILE=out.metrics ASO=coordinate VALIDATION_STRINGENCY=LENIENT
```

It works well but when the `input.bam` file gets bigger, the speed is very slow! I found that the picard MarkDuplicates doesn't support multiple threads. So, is there anyway to speedup picard? Another way, is there any better software to do the same as picard MarkDuplicates but with less time? I know [elprep](https://github.com/ExaScience/elprep) is another choice, but it needs very large memory!

Besides, I found that samtools can also remove duplicates, but according to my search, samtools can not remove the duplicates cross different chromosomes, so picard is more universe.

Any reply will be much appreciated!

Answer: > So, is there anyway to speedup picard?

1) picard stores data in memory until it needs to flush them on disk. The bigger the memory is (option MAX_RECORDS_IN_RAM) , the less you need I/O operation:

2) use sambamba-markdup http://lomereiter.github.io/sambamba/docs/sambamba-markdup.html (not tested)
Question: Before I report this as a bug can someone show me the proper way to add metadata using the `AddMetaData` function.

I first read a file that is new line delimited with the meta data labels I want using read table. The ftimecell_columns.txt file looks like this:

"12wks"
"12wks"
"12wks"
"12wks"
"12wks"
...

So I do this command to make the dataframe:

timecell_col <- read.table("~/Projects/FetalPancreas/timecell_columns.txt", header=FALSE, sep="\n")

and then this to AddMetaData:

AddMetaData(object = scfp, metadata = timecell_col, col.name = 'time.cell')

However, when I do:

head(x = scfp@meta.data)

All I see is the original meta data values:

orig.ident nCount_RNA nFeature_RNA
GSM2978830 scfetpan 382565 5080
GSM2978831 scfetpan 634726 4932
GSM2978832 scfetpan 565912 4501
GSM2978833 scfetpan 717152 5117
GSM2978834 scfetpan 207659 3508
GSM2978835 scfetpan 869148 5055

Am I reading the data in correctly?

Any help would be appreciated!

Very Respectfully,
Pratik

Answer: Is that your exact command line?

Not

scfp <- AddMetaData(object = scfp, metadata = timecell_col, col.name = 'time.cell')
Question: I am having some difficulty producing a script for parsing an alignment file in the following format generated from RepeatMasker. For example:

665 28.45 2.93 5.02 g5129s420 7350 7882 (1924) C MIR#SINE/MIR (1) 261 28 3

g5129s420 7350 ATCATAACAAACATTTAT--GGTGCCTCCTATGGAGCAGGGATTTTGCTT 7397
v v i i i v viv v i v v v
C MIR#SINE/MIR 261 ATAATAACCAACATTTATTGAGCGCTTACTATGTGCCAGGCACTGTTCTA 212

g5129s420 7398 AGGACTCTGAACTATAT---CTTACTT-GTCTTCATTAAAAACCTTATGA 7443
vi i iv i i i i i i v i
C MIR#SINE/MIR 211 AGCGCTTTACA-TGTATTAACTCATTTAATCCTCA-CAACAACCCTATGA 164

g5129s420 7444 AAAAGGTACTATTATTAACTGGGGXTGGGTTGTTTAACAGATAAGAAAGC 7787
iiv v i iii v i i i
C MIR#SINE/MIR 163 GGTAGGTACTATTATTATCC---------CCATTTTACAGATGAGGAAAC 123

g5129s420 7788 TTAAGAATTAGAGAGATAAATTATCTTGCTTAAGGTAACACAGTTAACAA 7837
v i v i i v v v ii v i ii
C MIR#SINE/MIR 122 TGAGGCA-CAGAGAGGTTAAGTAACTTGCCCAAGGTCACACAGCTAGTAA 74

g5129s420 7838 GCATTAG-GTCAAAGTTTGAACTCGGGCAGTCTGACTACAGAGCCC 7882
iivi i iiii i i i i v i
C MIR#SINE/MIR 73 GTGGCAGAGCCGGGATTCGAACCCAGGCAGTCTGGCTCCAGAGTCC 28

Transitions / transversions = 1.96 (45 / 23)
Gap_init rate = 0.03 (8 / 234), avg. gap size = 2.38 (19 / 8)

I would like to parse the file using BioPython such that I obtain the chromosome/scaffold name g5129s420, start + end (7350 7882), and the Transitions/transversions. Any ideas on how to write this script would be most welcome, as I am a complete novice to scripting.

Answer: There may be some BioPython modules or some other parser exists for this job. Here is my solution based on `Python 3.x` if you know `chromosome/scaffold name`

seqname = ['g5129s420','g5129s421'] #list sequence names here
dic, trans = {}, {}
for names in seqname:
dic[names] = []
trans[names] = []
with open("file.txt") as f:
for line in f:
line = line.rstrip().lstrip().split(' ')
line = [item for item in line if item is not '']
if len(line)>0:
#print(line)
if line[0] in seqname:
dic[line[0]].append(line[1])
dic[line[0]].append(line[-1])
seqname1 = line[0]
elif line[0] == 'Transitions':
trans[seqname1].append(line[4])
trans[seqname1].append(line[5].strip('('))
trans[seqname1].append(line[7].strip(')'))

for key, values in dic.items():
if len(dic[key])>0:
print(key,':', min(values), max(values), trans[key][0], trans[key][1], trans[key][2])

Output:

g5129s420 : 7350 7882 1.96 45 23
Question: I downloaded HapMap vcf file for one chromosome the file size is 11 Giga I am trying to use LiftoverVcf from Picard package, but it run out of memory. I am thinking to split the 11 giga file into smaller files each is one giga any tool can do that? or is there a better Idea?

already used ``` -Xmx200g -Djava.io.tmpdir=./mytmp ```

Answer: Using GATK:

- [RandomlySplitVariants][1]

>This tool takes a VCF file, randomly splits variants into different sets, and writes the results to separate files. By default the tool splits the input into two new sets, but it can be made to output more than two separate call sets.

[1]: https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_variantutils_RandomlySplitVariants.php
Question: I've been doing some work that involves characterizing potential chromatin states in a HCT116 cell model. I've successfully been able to run `ChromHMM` to identify chromatin states using a variety of histone markers and then overlap them with other datasets in order to double-check their annotation.

However, I am having a problem understanding some of the output files that `ChromHMM` automatically generates. Specifically `_emissions.txt` and `_*.bed`. I know there's a couple of people here that are really familiar with the software and could probably help me out.

I have already searched google, and read the ChromHMM manuscript, but neither provided answers.

Answer: The `_emissions.txt` are the values that go into the `_emissions.png` figures. Each row is a state, each column is an input data file ("mark" or histone mark in the terminology of ChromHMM). Darker blue indicates a higher likelihood of finding that mark in that state. These, combined with running OverlapEnrichment with biologically meaningful datasets, are critical for figuring out how to interpret the states.

The `segments.bed` file partitions the genome into contiguous segments, and the names of each feature in that file (E1, E2, etc) correspond to the states (1, 2, etc) in the `_emissions.png`.

A typical workflow is to figure out what to label each state. Then choose some colors and post-process the BED file with labels and names to get something more useful for downstream analysis.
Question: I have got a BED file (used in Plink and with .bim and .fam file) which contains Affy 500K data, I want to convert it back to PED file with map file. Is there anyone how to do this task? Thank you very much !

Answer: plink --bfile filename --recode --tab --out myfavpedfile

[plink](http://pngu.mgh.harvard.edu/~purcell/plink/dataman.shtml#recode)
Question: Hey ya'll,

Struggling here.

The dataframe was generated using a new package called InterCom ([https://github.com/saschajung/InterCom][1]), which basically determines ligand receptor pairing and also together with sending and receiving cells using scRNA-seq data. Incredibly useful if you ask me.

Here is a sample of my data:

"score" "significance" "lig.rec.pop" "ligand.receptor"
25507.7307692308 0.950689655172414 "undef_acinar" "FSTL1_BMPR1A"
113286.378762542 0.994068965517241 "mesenchymal_acinar" "FSTL1_BMPR1A"
56944.7025939177 0.989103448275862 "endothelial_acinar" "FSTL1_BMPR1A"
24467.0634475014 0.94551724137931 "mesenchymal_acinar" "BMP4_BMPR1A"
45149.7188075585 0.985034482758621 "mesenchymal_acinar" "BMP4_BMPR2"
20099.7619047619 0.916 "blood_acinar" "IL1B_IL1R1"
19786.641 0.913310344827586 "alpha_acinar" "VTN_ITGAV"
56609.9870823529 0.989034482758621 "mesenchymal_acinar" "FN1_ITGAV"

The maximum score value is 6611701.91, and the score minimum value is 36320.27,

Could someone help me get started with a good way to scale the data (if I even need to?, and generate a ComplexHeatmap? My main struggle is figuring out how to get it into the format to plug it into ComplexHeatmap.

I was able to do it through ggplot2, but I want to do it through ComplexHeatmap as I think it is easier to add bells and whistles. The only reason I was able to do it using ggplot2 was because I specified the what x, y and scores to use:

ggplot(hm, aes(x = ligand.receptor, y = lig.rec.pop, fill = score)) +
geom_tile() + theme(axis.text.x = element_text(angle = 90))

Ideally, I would love a heatmap/plot like this:

![enter image description here][2]

I realize that adding the time points (ie. week 7 and 9) might be a matter of adding meta data. but really if someone could just help me get started - I would be grateful.

Thank you in advance!

Image Citation: Gonçalves, C.A., Larsen, M., Jung, S. et al. A 3D system to model human pancreas development and its reference single-cell transcriptome atlas identify signaling pathways required for progenitor expansion. Nat Commun 12, 3144 (2021). https://doi.org/10.1038/s41467-021-23295-6

[1]: https://github.com/saschajung/InterCom
[2]: /media/images/1a0a4c91-f922-495b-8927-f5671604

Answer: Below chunk will convert your data to complexheatmap data format. U can then modify the heatmap.

Assuming your data is in **tdf** data frame.

lig = unique(tdf$lig.rec.pop)
rec = unique(tdf$ligand.receptor)

mat = matrix(0,nrow = length(lig), ncol = length(rec))
rownames(mat) = lig
colnames(mat) = rec

for(i in 1:nrow(tdf))
{
mat[tdf$lig.rec.pop[i], tdf$ligand.receptor[i]] = tdf$score[i]
}

Heatmap(mat)
Question: Hi All,
I need to split my fastq file that is composed of 30 chromosomes into 30 different files containing each the information from one chromosome. Technically, I need to split this kind of file into seperated one:

@chr1
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnatgctgggtgatctttagtcnnnnnnnnnn
nnnnnnnnnnnnnnnnatggggtcatgtacacacacacattggatannnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnatgctgggtgatctttagtcnnnnnnnnnn
...
@chr2
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnatgctgggtgatctttagtcnnnnnnnnnn
nnnnnnnnnnnnnnnnatggggtcatgtacacacacacattggatannnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnatgctgggtgatctttagtcnnnnnnnnnn
...

I tried to use :

awk '/^@chr1$/,/^+$/' consensus.fastq | perl -pe "s/@/>/ ; s/\+//" > chr1.fasta

but that gives me this :

@chr1
nnnnnagtnnnnnnnnnnnnnnnnnnnnnttgcnnnnnnnnnnnnnnnnnnnnnngcnnn
nnnntgaaannnnnnnnnnnntcnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn

Did somebody had already this kind of problems? Could somebody gives me some advices?
Thanks a lot!

Nath

Answer: If you have EMBOSS installed

seqretsplit -sformat fastq-sanger -osformat fastq file.fastq

will give you one Fastq record per file. The files will be named after the sequence ID, so chr1.fastq, chr2.fastq etc. Note that you should use `-sformat fastq-sanger`, `fastq-illumina` or `fastq-solexa`, depending on which encoding your file uses. See this question on Fastq format.

To get Fasta output (you don't explicitly ask for it, but your code implies it), simply change the `-osformat` argument.

seqretsplit -sformat fastq-sanger -osformat fasta file.fastq
Question: I grabbed the first couple of columns from my whole genome VCF file. It looks like this. Are there any tools either web-based, Python, or R that I could use to get SNP identifiers (e.g. `rs429358` or `rs7412`) for all of my SNPs that are in a particular database? I'm very new to working with VCF files and I want to figure out my blood-type. I'm very comfortable coding in both Python and R if there are any packages available for these languages. I would like to avoid, if possible, depositing my sequences in a 3rd party that would potentially use my information for their own gains but I am not opposed to having them as references in case I can't figure out any other options.

#CHROM	POS	ID	REF	ALT	QUAL
chrM	64	.	C	T	3070.00
chrM	73	.	A	G	3070.00
chrM	146	.	T	C	3070.00
chrM	153	.	A	G	3070.00
chrM	263	.	A	G	3070.00
chrM	310	.	T	C	3070.00
chrM	513	.	GCA	G	3070.00
chrM	663	.	A	G	3070.00
chrM	750	.	A	G	3070.00
chrM	1438	.	A	G	3070.00
chrM	1598	.	G	A	3070.00
chrM	1736	.	A	G	3070.00
chrM	1888	.	G	A	3070.00
chrM	2706	.	A	G	3070.00
chrM	3106	.	CN	C	3070.00
chrM	4248	.	T	C	3070.00
chrM	4769	.	A	G	3070.00
chrM	4824	.	A	G	3070.00
chrM	7028	.	C	T	3070.00
chrM	8027	.	G	A	3070.00
chrM	8794	.	C	T	3070.00
chrM	8860	.	A	G	3070.00
chrM	11719	.	G	A	3070.00

Answer: Make a BED file of SNPs for your genome of interest, e.g. `hg38` and dbSNP v147:

$ LC_ALL=C
$ wget -qO- http://hgdownload.cse.ucsc.edu/goldenpath/hg38/database/snp147.txt.gz \
| gunzip -c \
| awk -v OFS="\t" '{ print $2,$3,($3+1),$5 }' \
| sort-bed - \
> hg38.snp147.bed

Convert your VCF file to BED and map the ID field from the SNP file:

$ bedmap --echo --echo-map-id --delim '\t' <(vcf2bed --sort-tmpdir=$PWD < variants.vcf) hg38.snp147.bed > answer.bed

The `rs*` ID for the variant will be in the last column of `answer.bed`.
Question: I wonder what is the relationship between the depth and gene expression level in RNAseq?

For example, let's say I have RNAseq with 10M reads for one sample. If I use the same machine and generate 20M reads for the same sample, then except for the new detected genes, will the expression profile for genes expressed in both case be similar? In other words, if I plot the density VS log2(count) for them, will they have the similar shape?

Answer: Yes, if you use the same method to normalize the counts, you will get similar MA plots (in theory). At greater coverage, maybe you will see more dots, because there will be more genes detected.

The reads generated during the sequencing should be randomly distributed, except if the process runs into a problem. After normalization, you should have similar values for genes that were detected in both cases.

More coverage is useful to detect low expressed genes.

You will see significant differences in MA plot if you change the number of biological replicates. With more replicates, there will be more differentially expressed genes.
Question: Hi,

I would like to know how complex it is to analyze NGS Data. Is it possible to learn NGS data analysis from the online resources or should we learn under the guidance of an expert? How to get the core concepts of NGS data analysis? How to configure parameters while using open source tools? (Assembly, Alignment, statistics etc ). I have a masters degree in bioinformatics with unix, perl and basic core Java skills. Any advice is appreciated.

Answer: download bwa, samtools and a reference genome.

generate a random set of reads using samtools/misc/wgsim

index the genome

align the reads and generate a sam output.

describe each column of the sam

generate a vcf from the sam using samtools/mpileup

describe each column of the vcf

use ensembl/vep to predict the consequences of the variations.