text:
WGBS vs ONT, which one should I trust?
First of all let me apologize if this is not the place for this question. While it can be quite broad, I will try to make it as specific as I can.
Context
I am testing the correlation between Whole Genome Bisulfite Sequencing and Oxford Nanopore methylation calls for one individual (mammal, blood tissue) with replicates at all CpG sites. In the (very) best cases I find a close to 0.9 correlation. Yet, for the part that does not correlate I would like to know which technology is closer to the "biological truth".
My strategy
My plan would be to first find CpG sites/short CpG regions showing differential methylation states between WGBS and ONT. Then I would validate the methylation states of these sites using a third party technology such as pyromark.
What I am asking to the SE community

First, has anyone done this before ? I would love to hear about your experience

Do you see some critical aspects I should be watching out for ?

Last if you have comments about the strategy I propose here, you are welcome

I'd probably add EM-Seq (https://www.neb.com/products/e7120-nebnext-enzymatic-methyl-seq-kit) to the comparison, since it appears to become a valid replacement for WGBS based on enzymes rather than bisulfilte conversion.

@SebastianMüller Thank youfor the suggestion, I didn't know about EM-Seq so I will look at it.

Methylation can be cell-specific, which makes it difficult to evaluate accuracy on a bulk-cell level (even within the same tissue). How can you tell that the differences you're seeing are due to platform differences, or due to biological variation?
I find that adding more haystacks doesn't help much in working out the truth of a dataset. If you want to investigate biological truth, it would be better to create a biological system with a known methylation (or demethylation) pattern and test that.

thank you for your interest. Actually these datasets come from the same tissue sample. But you are right, working on a bulk-cell level makes it difficult to distinguish tech bias from biological variation. I think it can be difficult to produce control data that fits the biological variation you find in a bulk of wild cell, This is why I thought using a third party technology.

I would assess directionality and accuracy of prediction by 1) WGBS predicting ONT and then 2) ONT predicting WGBS.
Firstly, I would use deep learning (or machine learning) and train WGBS against ONT, parameterise and then test. Then conversely train ONT against WGBS, parameterise and then test. The approach with highest accuracy of prediction (using the 'accuracy' index) would be the approach assumed to be "closer to the biological truth". If both calculations produced comparable accuracies I would conclude what @gringer has stated that natural variation / natural heterogeneity is the predominant signal in the signal. This conclusion depends on the ability of deep learning to map the biological process.
It is an approach highly applicable to a deep learning estimate regardless and may not replace a truely controlled test, this will depend on what controls you have run, but could provide a valuable and easily obtainable insight in its absence.
If you did run a control sample, then you are truely in, because this would provide the training set and the WBGS and ONT provide the test data.

You asked me how it works
It doesn't replace good controls, but it is 'trendy workaround'.
Two steps

In deep learning (artificial neural networks, ANN, in this case you would use RNN such as LSTM) or machine learning (ML) e.g. random forest - random forests are MUCH easier to do, you supply a series of known values, thus CpG sites and their WGBS scores and use ANN (RNN approach) or ML (e.g. random forest). Both methods will look for pattens between the observed DNA data and the CpG occurence (training set). Usually this is 60% of your data set is used in training. A separate 20% of the data set is often used for parameterisation and a final 20% of the data set is used for testing. Thus how well can a ANN or ML predict a CpG site based on WGBS data. Likewise you would do the same for ONT. There is a caveat here which I discuss below.

The idea is if you give it a bit of DNA it will assign whether it is CpG (see caveat below).

However, instead of testing WGBS to WGBS you also test it WGBS to ONT and ONT to WGBS and the method with the highest 'accuracy' (positive hits, against true hits) is 'better method'. This is because you already know which is ONT positive and ask a WGBS trained algorithm whether it thinks its positive. The relationship between observed CpG and 'true (experimental) CpG is your accuracy score.

The caveat is that you need to supply the training set with sites that are known to be negative under WGBS or ONT. You need to think the best strategy for doing this. If some of those negatives were positve in the other method that would affect the outcome alot.
There is an issue about vectorisation this is often done using k-mers, and is a headache. Basically you have to make a bit of DNA into numbers and that in my opinion is the difficult part. In your case I ain't so sure its so difficult. If you can get that bit correct and biologically meaningful you are in business. My stuff is phylogeny based and I understand the relationship between biology (mutation) and numbers. In your case that needs thought, however someone will have thought about it and solved it.
As a personal statement I would do this via random forests, because ANN is difficult.
One such vectorisation method is here by Li et al (2017), however as they have used it, the method has flaws due to equal weighting of mutation frequencies. However, these flaws do not apply to your work - it would actually be a good method, they singly apply to tree building.

The final thing is that the more experimental data you can feed it the better. You don't have to complicated about it, just create a new column in the input section, e.g. along side teh vectorised DNA and the algorithm will figure out whether it helps. You can't feed it enough data, because it is doesn't help it wouldn't use it.

Thank you Michael, I am really interested by your approach but I think I don't understand what you mean by predicting WGBS (though I know the rôle of deep learning for ONT, I trained some models). I also have control samples (but the fully methylated one may not be trustable) .

@PaulEndymion I've modified the post above


text:
Cell to Cell communication detection tools for wide range of species
I am starting to get interested in protein interactome analysis from single cell/nucleus datasets.
I found three tools that might be interesting, which are CellPhoneDB, CellChat and ICellNet. However, these tools are limited to human or mouse cells and the developers' tips for using the tool for other species is to find the homologues or orthologues corresponding to the proteins in our species of interest.
What I understand is that these tools use their own database which has been curated using specific methods, which could mean that you would have to repeat the same steps for your species dataset to get a compatible database format and this could take a lot of time.
Also, I feel that using homologous proteins to find functional protein networks or protein-protein interactions would lead to many biases such as false positives due to the difference between human and other species proteins.
So I have a few questions here:

Has anyone already used this type of tool for species that the tool was not designed for? Did they find the results relevant?
Are there tools that allow the use of external databases (such as the STRING database)?
Perhaps you know of a similar tool that has been built specifically for Drosophila melanogaster (which is my species of interest).
Maybe you know of an alternative way/workaround to this problem ?

Thanks a lot!
Edit :

I found FlyPhoneDB which I am definitly going to test. However, I did not found any workaround for species with no dedicated tools. One reason is probably the time it would take to build such tool, though it would be highly valuable.
Closing this post because it too broad and generic.

This question is too broad and generic. Please limit your posts to a specific bioinformatics problem that you have encountered. If you want to ask multiple questions, then post multiple questions. If you want to have a general discussion about a topic, this is not the right place; ask your questions on a discussion forum like reddit.

Sorry, you are right. I thought about posting to reddit but I don't know much about the community, at least not the r/bioinformatics one.

Asking the question in the exact way you have here should be fine, with a flair of "discussion"; direct link to the post form here

I would recommend simply editing the question to focus on a single topic.You can ask multiple questions.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
samtools / bamUtil | Meaning of <wildcard> as Reference Name
I have been working with several tools on bam files recently and I am not sure how I should interpret some of the outputs.
samtools idxstats mapped.sorted.bam

Chr1    15096   0   0
Chr2    33397   0   0
Chr3    43888   41  0
Chr4    20819   1   0
*   0   0   34

Documentation says third and fourth column give "mapped reads" and "unmapped reads". So I interpreted the last line (*) as the number of unmapped reads as this is the only line with a non 0 value at this field.
Then I used
bam splitChromosome --in mapped.sorted.bam --out test/aln.

Reference Name: chr6 has 6 records
Reference Name: Chr7 has 1 records
Reference Name: Chr3 has 41 records
Reference Name: Chr4 has 1 records
Reference Name: * has 34 records

And looking at the bams outputted (one per reference sequence) in test/ I find an aln.UnknownChrom.bam that I am not sure where it comes from. I actually found 97 unique reads ID in it, where I was waiting for 34 or less.
So my first question is : Does "*" actually means "all reads that didn't mapp" and this information is stored in way I didn't know about in BAM files ?
And my second question would be : Anybody knows where this aln.UnknownChrom.bam comes from when using splitChromosome ? It is a source of error in a pipeline I am building and I'd like to get rid of it, if these are not actual mapping informations.

As per the SAM/BAM format specifications:

"RNAME: Reference sequence NAME of the alignment. If @SQheader lines are present, RNAME(if not ‘*’) must be present in one of the SQ-SN tag. An unmapped segment without coordinate has a ‘*’ at this field.  However, an unmapped segment may also have an ordinary coordinate such that it can be placed at a desired position after sorting.  If RNAME is ‘*’, no assumptions can be made about POS and CIGAR." -- http://samtools.github.io/hts-specs/SAMv1.pdf, emphasis mine

So, yes, these are the unmapped reads.
I haven't used this splitChromosome utility before myself, but as it splits the BAM file into one BAM file per chromosome, I suspect it simply also creates a BAM file for the unknown (*) chromosome. Why the numbers don't match I can't tell.
How to deal with that file in your pipeline will depend on what the pipeline is meant to do. But if you aren't interested in these reads, I suppose you could simply remove it before continuing to the next step.

Nice ! I'll maybe ask directly on bamUtil's github for my second question.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Where can I get ensembl Medaka genome for RNAseq
I am trying to map a RNA-seq dataset using ensembl genome for medaka fish.
From here http://uswest.ensembl.org/Oryzias_latipes/Info/Index
when I click on a Download DNA sequence (FASTA) I am forwarded here ftp://ftp.ensembl.org/pub/release-98/fasta/oryzias_latipes/dna/
where I have a huge number of .fa.gz files and I am just looking for a single file, that I can use.
While when I am downloading medaka genome from NCBI (and using this technique) http://www.metagenomics.wiki/tools/fastq/ncbi-ftp-genome-download I have my ncbi genome in one file.
How can I get an ensembl genome in one file?

Download the ones without a chromosome number in the filename. Those files are the whole genome.

deleted my answer, this is much better

Thank you Emily, but is that link located on a same page?? I mean the one " without a chromosome number in"? And you said oneS - there are multiple medaka genomes?

Yes, they're all listed in the same folder. Check the README to see what you want in terms of toplevel, repeat masking etc, but they're all there.

Thank you, but still... what I see is for example
Oryzias_latipes.ASM223467v1.dna_rm.primary_assembly.1.fa.gz Oryzias_latipes.ASM223467v1.dna_sm.primary_assembly.1.fa.gz``` but I cannot find something like ```Oryzias_latipes.ASM223467v1.fa.gz``` Sorry I am really new to bioinformatics and lots of things are new to me while they are somewhat obvious for the experts.. Thank you!

rm and sm refer to the repeat masking. rm is full repeat masking, all the repeats are replaced by Ns. sm is soft-masking, the repeats are in lower case. The other file has no repeat masking.


text:
Assign score to PyRanges from another PyRanges object
I have two PyRanges objects. PR1 is a set of regions of interest across a chromosome, PR2 tiles the genome and has a score associated with each range.
I need a way to get for each range in PR1 its score based on its overlap with PR2. If it is contained within just one range in PR2, it inherits the score of that range. If it overlaps multiple ranges on PR2, the score should be a weighted average of those ranges based on the size of the overlap.

Applying this to the image above, it would be:
score(A) = score(X)
score(B) = (OverlapBX * score(X) + OverlapBY * score(Y)) / sizeB

I've been browsing the PyRanges docs and I see a lot of functions that do almost that, but most of them don't take into account a score.
I had done this before in R using GRanges, where I could iterate over the overlap object to piece together scores, but I can't find a direct port of R's findOverlaps() in PyRanges.

The inputs I'm using for this are a list of mutation positions and a map of recombination rates over the chromosome they're on. So if two mutations are in a section with a constant recombination rate, the section between them inherits that recombination rate, but if the map says there's different rates between them, I need to take that into account.
Here's a pastebin with the data I'm working with: https://pastebin.com/sBYsjRQB

The Unfun Cat's first suggestion:
dist_to_pos = distances.join(positions, report_overlap=True).drop(like="_b")
df = dist_to_pos.df
g = df.groupby("Id2")
avg_score = g.cMperMb.sum() / g.Overlap.first()
avg_score.name = "AvgScore"
j_with_avg_score = pr.PyRanges(df.set_index("Id2").join(avg_score))

The results this give are... curious. Here's what it gives when applied to my data:
j_with_avg_score

Start   cMperMb     cMatstart   End     Chromosome  Id1     Overlap     AvgScore
0   2916415     1.223681    3.272236    2917484     chr1    0   1   12.236813
1   2916415     1.223681    3.272236    2917484     chr1    1   24  12.236813
2   2916415     1.223681    3.272236    2917484     chr1    2   5   12.236813
3   2916415     1.223681    3.272236    2917484     chr1    3   74  12.236813
4   2916415     1.223681    3.272236    2917484     chr1    4   110     12.236813
...     ...     ...     ...     ...     ...     ...     ...     ...
947     2960324     0.269208    3.291457    2960462     chr1    913     138     0.001951
948     2960462     0.352644    3.291494    2960790     chr1    913     77  0.018319
949     2960462     0.352644    3.291494    2960790     chr1    914     6   0.018319
950     2960462     0.352644    3.291494    2960790     chr1    915     15  0.018319
951     2960462     0.352644    3.291494    2960790     chr1    916     71  0.018319

952 rows × 8 columns
## input data contains 917 sections --> mismatch!

res_dic = {}
for x,y in zip(j_with_avg_score.Start, j_with_avg_score.AvgScore):
    res_dic[x] = y
res_dic

{2916415: 12.236813262999997,
 2917484: 0.02265139573888889,
 2917538: 0.06125764894,
 2917608: 0.10651243958260867,
 2918281: 0.6792060563708334,
 2918836: 1.1527867390640016,
 2921187: 4.176119918333334,
 2921274: 1.2543028034,
 2921300: 0.9058800047555555,
 2921976: 3.26162101122,
 2922306: 0.4522015679581395,
 2923593: 0.7409518310923077,
 2924172: 8.386433179800006,
 2925614: 0.31748233525000014,
 2926340: 0.17712928483599996,
 2926675: 1.4567026153333333,
 2926730: 5.063891720272731,
 2927172: 2.42991683355,
 2927410: 0.9370814839,
 2927411: 0.1303501307571428,
 2928125: 0.20534225654799998,
 2928557: 0.19228154158888885,
 2928849: 1.5796778519999994,
 2929269: 1.9946231915199983,
 2931597: 0.15826380113170732,
 2933043: 0.010025524,
 2934387: 0.05559167711694925,
 2938888: 0.0023231045055555556,
 2940194: 0.01882597932,
 2941104: 0.1670264778,
 2941694: 0.0197562495886076,
 2942700: 0.3931709166833331,
 2946526: 1.6788888500999994,
 2947460: 0.19597836537931007,
 2951262: 5.532616244500001,
 2956899: 0.44889279189999975,
 2960324: 0.0019507828376811595,
 2960462: 0.018319161953246753}

Something's getting chewed up here, due to the number of unique sections in the output, the duplicates, and the start and end positions not matching the input data.
However, the first value calculated is exactly the same as what I had gotten in my R script (with a few magnitude adjustments): 2916415: 12.236813262999997 vs 2917178: 0.00012236813263
So the calculations happening are most likely correct, I think it's mainly the start and end positions that are getting messed up.

I think join will get you halfway there? It is meant to be the replacement for findOverlaps

@TheUnfunCat I had started to fiddle with joins but then tried to use cluster and merge afterwards, but I got some numpy errors from those. I'll log them on the github.

@TheUnfunCat tried to replicate the error but it works now. I can't remember the details but it was something about numpy not recognising the variable type 'long'

Have you tried updating PyRanges?
Anyways, you can do gr.df to get the dataframe afterward your join and then do your calculations in pandas.

I have some inspiration below. It computes a weighted score based on overlap.
I need to know what exactly you want the overlap BX and BY to be. Could you post some example data and expected output?
# initialize test data
import numpy as np
np.random.seed(0)

import pyranges as pr
a1 = pr.data.aorta().drop()
a1.Id1 = list(range(len(a1)))
a2 = pr.data.aorta2().drop()
a2.Id2 = list(range(len(a2)))
a1.Score1 = np.random.random(len(a1)) * 100
a2.Score2 = np.random.random(len(a2)) * 100

# compute weighted average based on overlaps
j = a1.join(a2, report_overlap=True).drop(like="_b")

df = j.df

g = df.groupby("Id2")
avg_score = g.Score2.sum() / g.Overlap.first()
avg_score.name = "AvgScore"
j_with_avg_score = pr.PyRanges(df.set_index("Id2").join(avg_score))

Ideally I would like to combine BX and BY to B(XY), and get the average score during that process. I'll update the question with the files I'm using.

I tried the method, but it didn't work great. I'll update the post with the info from that


text:
Importing local files to a GATK docker image
I have pulled GATK docker on my computer but when running a GATK command, can not locate my input file
My OS version
cat /etc/os-release
NAME="Ubuntu"
VERSION="18.04.6 LTS (Bionic Beaver)"
ID=ubuntu
ID_LIKE=debian
PRETTY_NAME="Ubuntu 18.04.6 LTS"
VERSION_ID="18.04"
HOME_URL="https://www.ubuntu.com/"
SUPPORT_URL="https://help.ubuntu.com/"
BUG_REPORT_URL="https://bugs.launchpad.net/ubuntu/"
PRIVACY_POLICY_URL="https://www.ubuntu.com/legal/terms-and-policies/privacy-policy"
VERSION_CODENAME=bionic
UBUNTU_CODENAME=bionic

Files in my working directory
fi1d18@RBGO-Server2:/data/Continuum/WES/vcf$ ls
GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz
NG-27280_CLTSS_LTS_001A_lib506241_7636_2_MarkedDup.bam
NG-27280_CLTSS_LTS_001A_lib506241_7636_2_MarkedDup.bam.bai
NG-27280_CLTSS_LTS_001A_lib506241_7636_2_MarkedDuplicates.txt
NG-27280_CLTSS_LTS_001B_lib506242_7915_2_MarkedDup.bam
NG-27280_CLTSS_LTS_001B_lib506242_7915_2_MarkedDup.bam.bai
NG-27280_CLTSS_LTS_001B_lib506242_7915_2_MarkedDuplicates.txt
NG-27280_CLTSS_LTS_002B_lib488377_7494_1_MarkedDup.bam
NG-27280_CLTSS_LTS_002B_lib488377_7494_1_MarkedDup.bam.bai
NG-27280_CLTSS_LTS_002B_lib488377_7494_1_MarkedDuplicates.txt
fi1d18@RBGO-Server2:/data/Continuum/WES/vcf$

I then pulled GATK docker
fi1d18@RBGO-Server2:/data/Continuum/WES/vcf$ docker run -v ~/vcf:/gatk/vcf -it broadinstitute/gatk:latest
(gatk) root@8b2af4d7aadf:/gatk#

(gatk) root@8b2af4d7aadf:/gatk# ls
GATKConfig.EXAMPLE.properties  gatk-package-4.4.0.0-local.jar  gatkPythonPackageArchive.zip  scripts
README.md              gatk-package-4.4.0.0-spark.jar  gatkcondaenv.yml          vcf
gatk                   gatk-spark.jar              gatkenv.rc
gatk-completion.sh         gatk.jar                run_unit_tests.sh
(gatk) root@8b2af4d7aadf:/gatk#

As you see, folder seems to be empty because I can not connect my local working directory to container
(gatk) root@6d3d3ff69e21:/gatk# cd vcf
(gatk) root@6d3d3ff69e21:/gatk/vcf# ls
(gatk) root@6d3d3ff69e21:/gatk/vcf#

I have tried all of these but did not work
fi1d18@RBGO-Server2:/data/Continuum/WES/vcf$ docker run -v /data/Continuum/WES/vc:/gatk/ -it broadinstitute/gatk:latest

fi1d18@RBGO-Server2:/data/Continuum/WES/vcf$ docker run -v ~/data/Continuum/WES/vc:/gatk/ -it broadinstitute/gatk:latest

What does -v ~/vcf:/gatk/vcf do? Why are you using ~/vcf there? ~ is your $HOME directory, is /data/Continuum/WES/ your home directory? Does it work if you try docker run -v /data/Continuum/WES/:/gatk/vcf  -it broadinstitute/gatk:late?

Yes you solved that, thanks a lot

@terdon answered that
docker run -v /data/Continuum/WES/:/gatk/vcf -it broadinstitute/gatk:latest

Glad it's solved. Unless @terdon posts this as an answer, could you kindly mark this as "resolved" please (this helps with site stats). The "resolved" flag can always be shifted between posts, i.e. if terdon answers later.

I do not know how to do that sorry

Just "accept" your answer (see https://bioinformatics.stackexchange.com/help/someone-answers) by clicking on the checkmark on the left.


text:
A bash script for running on a bunch of bam files
I have some bam files in this directory
/data/Continuum/WES/results/

I want to run GATK mutation calling over bam files
I googled and I realised for one function, I can do this
cd /data/Continuum/WES/vcf/

for file in *.bam ; do ./gatk CollectSequencingArtifactMetrics -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup --FILE_EXTENSION .txt -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta done;

I have five functions to run but I do not know how to include all in the script
These are my functions
./gatk CollectSequencingArtifactMetrics -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup --FILE_EXTENSION .txt -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta

./gatk GetPileupSummaries -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.table -V af-only-gnomad.hg38.vcf.gz -L resources_broad_hg38_v0_wgs_calling_regions.hg38.interval_list -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta

./gatk CalculateContamination -I NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.table -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.contamination.table

./gatk GetSampleName -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.sample_name

./gatk Mutect2 -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.vcf -tumor NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.sample_name --af-of-alleles-not-in-resource 2.5e-06 --germline-resource af-only-gnomad.hg38.vcf.gz -pon 1000g_pon.hg38.vcf.gz

I thought to use | to include five functions but I do not know how to change input and output for each iteration
I tried this for the first function,
    FILES="/data/Continuum/WES/testAlignmentBROADGenome2/results/.MarkedDup.bam*"
for f in $FILES
do
  ./gatk CollectSequencingArtifactMetrics -I $f -O $f --FILE_EXTENSION .txt -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta
done

but returned this error
SAMException: Cannot read non-existent file: file:///data/Continuum/WES/testAlignmentBROADGenome2/results/.MarkedDup.bam*

Anyone, please help me in develop this?

There are a couple ways to approach this, one is to create output files based on the BAM input file name - the other approach is to make a new folder (based on the input BAM filename containing result files with the same name for each input. I generally use the first method, but either work.
Here is a script that may run the 5 GATK functions. It calculates a prefix to be used for derivative output files BASE, and should run your pipeline on all the BAM files within the INPUT_FOLDER.
#!/bin/bash

set -e

INPUT_FOLDER="/data/Continuum/WES/results/"

for BAM_FILE in "$INPUT_FOLDER/"*.bam
do
        BASE=$(basename "$BAM_FILE" .bam) # NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup

        # function 1
        ./gatk CollectSequencingArtifactMetrics \
                -I "$BAM_FILE" \
                -O "$BASE" \
                --FILE_EXTENSION .txt \
                -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta

        # function 2
        ./gatk GetPileupSummaries \
                -I "$BAM_FILE" \
                -O "$BASE.targeted_sequencing.table" \
                -V af-only-gnomad.hg38.vcf.gz \
                -L resources_broad_hg38_v0_wgs_calling_regions.hg38.interval_list \
                -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta

        # function 3
        ./gatk CalculateContamination \
                -I "$BASE.targeted_sequencing.table" \
                -O "$BASE.targeted_sequencing.contamination.table" \

        # function 4
        ./gatk GetSampleName \
                -I "$BAM_FILE" \
                -O "$BASE.targeted_sequencing.sample_name"

        # function 5
        ./gatk Mutect2 \
                -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta \
                -I "$BAM_FILE" \
                -O "$BASE".vcf \
                -tumor "$BASE.targeted_sequencing.sample_name" \
                --af-of-alleles-not-in-resource 2.5e-06 \
                --germline-resource af-only-gnomad.hg38.vcf.gz \
                -pon 1000g_pon.hg38.vcf.gz
done

It may need some adjustments - I did not test this, and used the GATK incantations you provided (thank you for including all of your steps!)
It includes the addition of the command set -e, which depending on GATKs return status may stop the script if any command fails - otherwise the script will happily run through a failed command.
You also may need to adjust the location of your output, and the directory the script is in upon execution, as you are running GATK from the current directory (./gatk)
In order to make these long lines readable, I use the method of placing a \ at the end of a line - this allows the command to wrap to the next line. Be sure that the \ character is at the end of your line for this to function properly.

Wow, cool, thanks.

OMG Thanks a million the script is running

Sorry @Scot, as the Mutect2 part of the script takes 15 hours in average for each sample, does this mean that for function 5, script will wait for the first sample to being completed by Mutect2 then goes for the second sample or Mutect2 will run of all of 18 samples in parallel? In the first case, the script will take 18 days?

Shell scripts are run line by line. If a line takes 15 hours, the interpreter will not move to the next line before that line is done. Loops are merely instructions to run a set of lines multiple times in series, so there's no parallel action happening here. You should ideally write a snakemake script with a rule for each GATK step then run a bunch of snakemake processes in parallel using whatever parallelization means are available to you.

Hi @Angel - yes, the current version of this script runs sequentially, and is not parallelized. If function5 Mutect2 takes a long time, then it may be worthwhile to pull it out of the loop, and run it upon completion of all other steps, and addition, explore if it can be run in parallel.

You really don't want to use CAPS for shell variable names. Since, by convention, global environment variables are capitalized, it is bad practice to also use caps for your own, local variables because this can lead to naming collisions and hard to find bugs.


text:
Make a new column in R or bash
I have antibody data the first column is heavy chain sequence and the second column is the light chain
I want to create a third column by combining two columns like
If CARGRGIAVAGTTGFDYW is a heavy chain and CQQYNSYSLTF a light chain, in third column I want CARGRGIAVAGTTGFDYW_CQQYNSYSLTF
if a heavy or light chain is absent, I want something like NA_CQQYNSYSLTF or CARGRGIAVAGTTGFDYW_NA depends of absent part
CDR3.light.aa   CDR3.heavy.aa
-              CAGDLSVGGVLAGGVPHLRHFDPW
-              CPETPRYCSGGFCYGEFDNW

CSSYAGSNTFGLF   -
CSSYAGSNTGLF    -

So far I have done this
a$a <- paste(a$CDR3.heavy.aa, a$CDR3.light.aa, sep= "_")

But I noticed a lot of mess and manual works
Any help please?

What did you try so far?

Thank you, I have added my code in the post

Your approach was correct, you just have to transform your dataframe first to replace the character "-" by NA :
library(stringr)
library(dplyr)
df=data.frame(CDR3.light.aa=c("-","-","CSSYAGSNTFGLF",'CSSYAGSNTGLF'),
              CDR3.heavy.aa=c("CAGDLSVGGVLAGGVPHLRHFDPW","-","CPETPRYCSGGFCYGEFDNW","-"))

df=df %>%
  mutate_all(funs(str_replace(., "-", "NA")))

df$a <- paste(df$CDR3.heavy.aa, df$CDR3.light.aa, sep= "_")

 df
  CDR3.light.aa            CDR3.heavy.aa                                  a
1            NA CAGDLSVGGVLAGGVPHLRHFDPW        CAGDLSVGGVLAGGVPHLRHFDPW_NA
2            NA                       NA                              NA_NA
3 CSSYAGSNTFGLF     CPETPRYCSGGFCYGEFDNW CPETPRYCSGGFCYGEFDNW_CSSYAGSNTFGLF
4  CSSYAGSNTGLF                       NA                    NA_CSSYAGSNTGLF


text:
Creating a scatterplot by a json data in js
I have a json file of protein data of mouse and human like
{
    "name": "P04202",
    "full_name": "Transforming growth factor beta-1 proprotein",
    "symbol": "TGFB1_MOUSE",
    "length": 390,
    "species": "mouse",
    "function": "Transforming growth factor beta-1 proprotein: Precursor of the Latency-associated peptide (LAP) and Transforming growth factor beta-1 (TGF-beta-1) chains, which constitute the regulatory and active subunit of TGF-beta-1, respectively.",
    "nr_modifications": 2821,
    "residues": [
        {
            "pos": 1,
            "residue": "M",
            "modifications": []
        },
        {
            "pos": 2,
            "residue": "P",
            "modifications": []
        },
        ...
        {
            "pos": 60,
            "residue": "A",
            "modifications": [
                {
                    "mod": "[127] Fluoro",
                    "type": "Chemical derivative",
                    "prob": 1
                },
                {
                    "mod": "[53] HNE",
                    "type": "Post-translational",
                    "prob": 0.25
                }
            ]
        },
        ...
        ]
}

How do I create a scatter plot [![enter image description here][1]][1]
Thanks for any idea
So far I gathered a code to fetch the data from an online json
export const load = ({ fetch }) => {
  const fetchProteins = async () => {
    const res = await fetch('https:address/assets/ptm.json')
    const data = await res.json()
    return data
  }

  return {
    Proteins: fetchProteins()
  }
}

This code generates a scatter plot but I will need to hover around the dots
script>
    import { extent } from 'd3-array';
    import { scaleLinear } from 'd3-scale';
    export let data = [];

    const margin = 30;

    $: xScale = scaleLinear()
        .domain(extent(data.proteins.map((d) => d.length)))
        .range([0, 800 - (margin * 2)]);

    $: yScale = scaleLinear()
        .domain(extent(data.proteins.map((d) => d.nr_modifications)))
        .range([0, 400 - (margin * 2)]);
</script>

<svg width="800" height="400">
    <g transform="translate({margin}, {margin})">
        {#each data.proteins as datapoint}
            <circle
                cx={xScale(datapoint.length)}
                cy={yScale(datapoint.nr_modifications)}
                class:human={datapoint.species === 'human'}
                r="20"
            />
        {/each}
        <circle
                cx={10}
                cy={10}
                class="human"
                r="10"
            />
            <circle
            cx={10}
            cy={40}
            class="mouse"
            r="10"
        />

    </g>
</svg>

What did you try so far?

I gathered a code to fetch the data from an online json

"fetch data from online JSON" is basically downloading a file and maybe subsetting it. What else have you tried in actually generating the plot?

Yesterday I accomplished to plot a scatter chart but I will need to deploy hover around the dots

It seems like you are already on the right track. The code you provided fetches the protein data from the JSON file and generates a scatter plot with D3.js based on two properties: length and nr_modifications.
To add an interactive hover event, you can add a tooltip to your scatterplot. The tooltip will display additional information when you hover over each dot.
Below is an example of how you can add a tooltip:
First, add a
 element in your HTML that will be used as the tooltip:

Protein:

Modifications:

Length:

Then, in your script, add the code to handle the mouseover and mouseout events for each circle:
{#each data.proteins as datapoint}
    <circle
        cx={xScale(datapoint.length)}
        cy={yScale(datapoint.nr_modifications)}
        class:human={datapoint.species === 'human'}
        r="20"
        on:mouseover="{(e) => handleMouseover(e, datapoint)}"
        on:mouseout="{(e) => handleMouseout(e)}"
    />
{/each}

Define the handleMouseover and handleMouseout functions in your script. These functions will show and hide the tooltip respectively. For the handleMouseover function, it also updates the tooltip content with the protein data:

In standard HTML or JavaScript, the {#each} syntax and usage of {} for expressions and variable binding won't work as they are specific to Svelte. Same applies to the usage of on:mouseover and on:mouseout. In standard HTML or JavaScript, we use onmouseover and onmouseout.
Here are some things which I haven't checked, but you might want to consider looking at:

It looks like you are trying to load JSON data from an external
source, but the function fetchProteins() is never called and its
result (the promise returned) is not being awaited. The way it's
structured now, Proteins in the returned object will hold a Promise,
not the actual data.
The code assumes that the returned JSON data is an array stored in
data.proteins, but the provided JSON example suggests that the data
structure might be different. It doesn't show any 'proteins' key in
the JSON. The code needs to align with the actual structure of the
returned JSON data. data.proteins implies that data is an object and proteins is a property of that object containing an array.
In your SVG transform attribute, you are using braces {} though the correct syntax for a
translate transformation in SVG is translate(x,y).

You might want to try looking at the following changes:
export const load = ({ fetch }) => {
  const fetchProteins = async () => {
    const res = await fetch('https:address/assets/ptm.json')
    const data = await res.json()
    return data
  }

  return {
    Proteins: fetchProteins
  }
}

  $: xScale = scaleLinear()
        .domain(extent(data.map((d) => d.length)))
        .range([0, 800 - (margin * 2)]);

    $: yScale = scaleLinear()
        .domain(extent(data.map((d) => d.nr_modifications)))
        .range([0, 400 - (margin * 2)]);

<g transform="translate({margin}px, {margin}px)">
        {#each data as datapoint}
            <circle
                cx={xScale(datapoint.length)}
                cy={yScale(datapoint.nr_modifications)}
                class="{datapoint.species === 'human' ? 'human' : ''}"
                r="20"
            />
        {/each}

</g>


text:
Variant allele frequency versus population frequency
In gnomAD related article
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9160216/#:~:text=Allele%20frequency%20is%20calculated%20by,can%20differ%20substantially%20between%20positions
I road this section

Allele frequency is calculated by dividing allele count by allele
number, hence, allele frequency represents the frequency of
confidently sequenced haplotypes that carry the allele in question
(because coverage and sequencing quality varies across the genome, the
allele number can differ substantially between positions). Allele
frequency is not equivalent to the percentage of individuals that
carry the allele, but is a suitable value for expressing the frequency
of a variant in the general population. The number of individuals
carrying a variant will depend on the number of heterozygous and
homozygous individuals but can be calculated from the data provided in
the variant table.

Does this mean, for getting the frequency of each allele in gnomAD, one could use Allele frequency (AF) column or must calculate population frequency?

What VCF fields are you comparing? I see the gnomAD has the overall AF field and various AF fields calculated by sample subset.

How did you subset gnomAD VCF to non-finnish samples only?

You might find gnomAD's non-cancer AF fields useful. I'll try to decipher the formula. AN/2 is the number of samples (assuming all samples are diploid). AC - num_homo gives you 2*(num(hets)). Dividing this by number of samples gives you frequency of samples with at least one ALT allele.

gnomAD provides a global Allele Frequency as well as AF values for most subsets. gnomAD does not provide individual genotypes (as far as I know), so calculating custom AF for non-precomputed subsets might not be possible.
Since you're interested in the NFE (non-Finnish European) sub-population, the nfe_*_AF fields might be of interest to you.

As far as calculating the fraction of non-WT samples from num_homozygous_samples, AC, and AN goes, your formula seems spot on. Here's my deconstruction of the formula:
Your formula is: (AC - num_homozygous_samples)/(AN / 2))
AN/2 is the number of samples (assuming all samples are diploid). AC - num_homo gives you 2 * num_hets. Dividing this by number of samples gives you frequency of samples with at least one ALT allele.
For example, let's take a dataset of 200 diploid samples where 100 are WT, 75 are het and 25 are hom-alt for a biallelic variant.
Since there are 200 diploid samples, AN is 400.
AC for the ALT allele is (75 * 1) + (25 * 2) = 125.
Since 100 are WT, it is known that fraction of samples with mutation = 1 - (100/200) = 0.5. This number can be arrived using your formula:
(AC - num_homozygous_samples) / (AN/2)
 = ((125 - 25) / (400/2))
 = 100 / 200
 = 0.5


text:
Population Allele Frequency
In UK biobank research environment
In GENOMICS tab of UK BRAP, if I type 7:140476600-140480852 coordinate in search bar and the same time after filtering my cohort for this coordinate, in ALLELE TABLE I see columns like:
Cohort AF Cohort Allele Count Population AF Population Allele Count GnomAD AF

I Googled a lot but no clue where population allele count and frequency come from
Does anyone have any experience with this?

I don't know specifically about the UKBB RAP, however, Population Allele Count is the count of that specific allele genotyped in the population in question. In this case, it most likely refers to the population that your cohort is part of. The population allele frequency is the same thing but noted in terms of frequency, i.e. divided by 2 times the number of genotyped individuals in that population.


text:
Printing the results of for loop into a txt file
I have a for loop which prints available vcf files in the path for a list of patient IDs (using find dx data function) into my zsh terminal:
for i in $(cat patient_id.txt); do
 dx find data --property patient_id=$i --path "vcf/"
done

How I can have a file contain the vcf files that the loop shows in terminal? I mean by using this command in the terminal, I am seeing the name and path of my vcf files for a given participant list and working directory; Instead of seeing the file name in the terminal, I want to have that file name in a txt file.
Thanks

Please [edit] your question and explain what happens if you simply add > file after the done. Is output redirected to file? Is it still shown in your terminal? What is the output? Are you sure files are found and you aren't just getting error messages?

This command worked just by adding  > file after done
for i in $(cat patient_id.txt); do
dx find data --property eid=$i --path 'vcf/'
done > vcf_list_file.txt

Things that write to standard output on the terminal can be wrapped up in brackets and redirected to a file. Assuming this is what you want to do, I added brackets around the existing code, then used > to redirect the output to a file.
(for i in $(cat patient_id.txt); do
 dx find data --property patient_id=$i --path 'vcf/'
done) > vcf_files.txt

Note that the brackets go around the entire code, including the for bit.

@zizogolu the brackets go around the entire code, including the for bit

@Zizogolu can you please type out what you put into the terminal, and the error messages that were produced?

There is no need to make this run in a subshell (that's what the parentheses do). All you need is the redirect in the end. for i in ...; do ...; done > file works and redirects all stdout of the for loop. Compare for i in a b; do echo "$i"; done > file and the same thing with parentheses, it makes no difference other than running the whole thing in a subshell.

Thanks a million @terdon , your comment solved my issue

That is... odd, @Zizogolu. If it works with for ...; do ... done > file then it would also work with Gringer's (for ...; do ... done) > file, the parentheses are just irrelevant, they wouldn't stop the redirection from working. In any case, can you please post what worked for you as an answer so the question can be marked as answered?

I suspect the issue was in the path command-line argument, which seems to have changed since I copied the code. I've updated that now.


text:
Which data is being used for violin plot?
Sorry I just got totally confused conceptually.
If this is my raw count data
Likely, Seurat divides each value by sum of the column afterward times by 10000.
which gives so
                   s1.1        s1.2 s1.3       s1.4 s1.5
DDB_G0267178 0.00000000 0.009263254    0 0.01286397    0
DDB_G0267180 0.00000000 0.000000000    0 0.00000000    0
DDB_G0267182 0.00000000 0.000000000    0 0.03810585    0
DDB_G0267184 0.00000000 0.000000000    0 0.00000000    0
DDB_G0267188 0.02640801 0.000000000    0 0.01286397    0

Literally Seurat claims that final step would be taking natural log of the above matrix as I tested with
library(seurat)
mat <- matrix(data = rbinom(n = 25, size = 5, prob = 0.2), nrow = 5)
mat

    [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    1    1    1
[2,]    1    0    2    2    3
[3,]    1    1    1    1    2
[4,]    1    1    0    1    0
[5,]    0    0    0    1    3
mat_norm <- LogNormalize(data = mat)
mat_norm
[1,] 7.824446 .        7.824446 7.419181 7.014015
[2,] 7.824446 .        8.517393 8.112028 8.112028
[3,] 7.824446 8.517393 7.824446 7.419181 7.706713
[4,] 7.824446 8.517393 .        7.419181 .
[5,] .        .        .        7.419181 8.112028

My confusion is: seurat@data dose not give me natural log transformed data rather returns only data divided by column sum and scalled by 10000.
So when plotting a violin plot if data are log normalised plot should not look as below picture
When I used
VlnPlot(object = y, features.plot = "DDB_G0277853", x.lab.rot = FALSE, y.log = TRUE)

plot changed totally. Could someone please help me in getting idea about seurat@data that is not really log transformed and by which data violin plot is being produced??? So if a gene has 10000 read counts, without log transformation violin plot should be so
This is my seurat object

The normalized and log-transformed values are used for the violin plot.
The argument y.log changes only the display of the data (scaling of the y axis).
Seurat has very good documentation. Section 7 in the FAQ explains what data is stored in the object:

How is data stored within the Seurat object? What is the difference between raw.data, data, and scale.data?

raw.data The raw data slot
(object@raw.data) represents the original expression matrix, input
when creating the Seurat object, and prior to any preprocessing by
Seurat.
data The data slot
(object@data) stores normalized and log-transformed single cell
expression. This maintains the relative abundance levels of all genes,
and contains only zeros or positive values. See ?NormalizeData for
more information. This data is used for visualizations, such as violin
and feature plots, most differential expression tests, finding
high-variance genes, and as input to ScaleData (see below).
scale.data
The scale.data slot (object@scale.data) represents a cell’s relative
expression of each gene, in comparison to all other cells. Therefore
this matrix contains both positive and negative values. See ?ScaleData
for more information.

Thank you, you are definitely alright. The problem is; when I am normalising values for some genes and compare them with output of seurat@data, some genes are ok, but some genes seems to be without log normalisation and all stored in seurat@data. I guess maybe after calling by 10000, only some especial range of expression values go through the log normalisation, for example seurat@data gives me 0.0026 for a gene in a cell, when I calculated  by pen, I noticed this value is just deciding the raw count on column sum multiply by 10000. Please correct me if I am wrong

In your first example, the last row does not match: DDB_G0267186 vs DDB_G0267188. Note that object@raw.data contains all cells, meanwhile object@data contains only the filtered cells. That might explain the difference.

@FereshTeh look for comparison at the formula Seurat uses: log1p(value/colSums[cell-idx] *scale_factor) (note the implicit +1 and logn in the formula)

Those two graphs don't look radically different.   The y axis is definitely not log scaled on the top one, and it looks like it is scaled in the second one (y.log = TRUE)

Sorry, what you think about logNormalize function in Seurat that does not return log  transformed data?

Sorry, for instance, for gene 1 in cell 1 I have 1040 RAW read counts. the sum of raw read counts in cell 1 for all genes is 373762. Based on the seurat documentation for normalisation, I should do so;  ln ( 1040/373762 * 10000) that gives me 3.6 that is right. So in top violin plot for this gene seurat has used the natural log transformed data as stated and in seconds violin plot I am just taking another log of y axis. But question is: why seurat@data does not return values as I expect and I called for one gene rather give me only scaled data by 10000. could you help me?


text:
Violinplot of gene expression
I have plotted the log normalized expression of two genes by violonplot for 4 clusters. I have links to my pictures and Seurat object too. My problem is this; in violin plot I can not see the mean or any centennial tendencies so that I don't know if two genes is expressing higher or lower in contrast to each other in each cluster.
violonplot vlnPlot(seuset,c(""))
Actually when I am thinking I see I need a violin or box plot showing  the  ratio of these genes to each other in each cluster and how nice to have mean also in the plot
Somebody knows how I can do that?

Please, remember to add the code you use to make it easier to provide the accurate advise to help you. But do you want to see the mean of the cluster or to see the differences of genes between clusters?

You want to (1) see the mean for each gene, and also to (2) calculate a ratio of expression levels of two genes, then compare it between clusters.
(1) First, notice that vlnPlot() is deprecated. Use VlnPlot().
You can try using the parameter do.sort=T:
VlnPlot(object=seuset, features.plot=c("DDB_G0267412", "DDB_G0277853"), do.sort=T)
Alternatively, you can return the ggplot2 object and then plot the means etc:
p <- VlnPlot(object=seuset, features.plot="DDB_G0277853", do.return=T)
p <- p + geom_boxplot(width=0.05)
p

(2) There are a few problems with calculating the ratio of gene expression levels. Firstly, what do you mean by gene expression level and how do you measure it? Surely you can have the read counts, but how do you interpret them? Is 30 counts a high level? Is it much more than 60 counts, or is it roughly the same? The same applies to the calculated ratios and the differences between them, even if we ignore amplification, gene length and other biases. And it is very hard to interpret ratios if the reference can also change. The problem is somewhat similar to RT-qPCR, where people use a set of reference genes whose expression has previously been shown to be invariant under the conditions.
In theory, you could use the raw counts (object@raw.data), the log + normalized counts (object@data), or the scaled counts (object@scale.data). The raw counts are biased by sequencing-depth, and the ratio of log or scaled values are not easily interpretable or intuitive. Therefore the library-size normalized (non-log) values seem to be the best.
Additionally, you could calculate the ratio of two genes either (a) for each cell (paired), or (b) for each group.
(a) is problematic, because of the zero values: you will have many NaN and Inf values, which cannot be removed without biasing the data.
Here are some functions for retrieving and plotting data from the object:
ident1 <- WhichCells(seuset, ident=1)
ident2 <- WhichCells(seuset, ident=2)

FetchData(seuset, vars.all='DDB_G0267412', cells.use=ident1, use.raw=T)

datatoplot1 <- seuset@raw.data['DDB_G0267412', ident1] / seuset@raw.data['DDB_G0277853', ident1]
datatoplot2 <- seuset@raw.data['DDB_G0267412', ident2] / seuset@raw.data['DDB_G0277853', ident2]

boxplot(datatoplot1, datatoplot2)

stripchart(datatoplot1,
            vertical = TRUE, method = "jitter",
            pch = 21, col = "dark green",
            add = TRUE, at=1)

stripchart(datatoplot2,
            vertical = TRUE, method = "jitter",
            pch = 21, col = "dark green",
            add = TRUE, at=2)

excuse me, with this command i have this picture in the link for my four clusters p <- VlnPlot(object=seuset, features.plot="DDB_G0277853", do.return=T)
p <- p + geom_boxplot(width=0.05)   https://drive.google.com/file/d/1r6eGQB225_jwtf7AWQo6Z2FKj4eVro42/view?usp=sharing, how i can have the plot for overall data side by side this plot?

i plotted that for all of cells but i don't know how to make a 5 violin together. https://drive.google.com/file/d/1MarsjXbTf0jg8e8e-1MTARNFQsOC9svw/view?usp=sharing   so i would be able to compare each cluster mean with the overall  mean, how i can have the plot for overall data side by side this plot? i plotted that for all of cells but i don't know how to make a 5 violin together


text:
Why DoHeatmap Does not show all genes in genes.use?
I am heatmaping a list of genes by DoHeatmap function in Seurat R package. I am sure I have 212 genes but heat map shows only a few of my genes
> DoHeatmap(
+        object = seurat,
+        genes.use = genes,
+        slim.col.label = TRUE,
+        remove.key = TRUE)
> dim(as.matrix(seurat@data))
[1] 12293   209

> length(genes)
[1] 212
> class(genes)
[1] "character"

> class(seurat)
[1] "seurat"
attr(,"package")
[1] "Seurat"
>

Please look at heat map, I was trying to plot 212 genes but you seeing only few genes are among my 212 genes
Please somebody save me from this confusion

This is probably a very naive question, but do more genes appear if you make the image larger? I mean, if you extend the image before exporting it, while it is being displayed in R?

Actually no, because before exporting from R looks the same :(

Yes, but what if you resize it before exporting? Does making the image bigger make more genes appear? I just want to rule out the obvious here.

Thank you, I just change the size of image before exposing to 2500*1991 but only rows became smaller and the order and number of genes are the same

Oh well. It was not very likely to be that simple.

https://github.com/satijalab/seurat/issues/287
In this post there is a answer for my question, however when I set scale to FALSE again the number of genes did not change but the shape of plot dramatically changed. So, maybe the issue is because of scaling that did not solve my problem


text:
Calculating L2 norm (Euclidean norm) to set of genes in time series
I have gene expression RNA-seq data of 9 time points (0 hour, 2 hour, ..., 16 hour). I need a plot on which x axis is time point and each data point at a certain time point represents the L2 norm (Euclidean norm) of the difference between the gene expression vector at that time point and the previous time point.
This is the link of my expression data
https://drive.google.com/file/d/1crJSqVdGujy1tFyeoV6E7e-xO9SXteJZ/view?usp=sharing

If you read your matrix into R and saved it as m:
diffMatrix = m[,2:9] - m[,1:8]  # The change per-timepoint
l2n = apply(diffMatrix, 2, function(x) sqrt(sum(x^2)))  # The L2-norm

R also has a norm() function, which can directly produce the L2-norm, but it turns out that this is slower than the simple method shown here. If you click through to the other question you can find somewhat faster implementations of this, but anything will suffice for small datasets like yours.

Sorry, How can I normalise these 8 points to zero time point?     12n= Named num [1:8] 313698 168553 360237 487963

That depends entirely on how you want the normalization done.

Sorry I mean how I can normalize these 9 differences you helped me to calulate  to expression at 0 h, which corresponds to cells growing in ground state conditions?? Because they are very big and I can not plot them

You wouldn't, it doesn't make sense to normalize a difference between two time points to a third one. You could normalize everything before subtraction though.

Right, actually the goal was showing global gene expression dynamics as measured by bulk RNA-seq. I was going to plot a data point at a certain time point represents the 2-norm of the difference between the gene expression vector at that time point and the previous time point. However these differences are too big to be plotted. Thank you

Try something like plotting log2(TPM) in violin or joy plots. If there are changes in the distribution of expressed genes over time then that should make it apparent.


text:
Adding metadata to Seurat object
I have a Seurat object of 8 patients. I want to add metadata to that so that I have origin of each cell. At the moment UMAP just shows a bunch of cells while I want to color clusters by sample
This is my metadata
> str(head(metadata))
'data.frame':   6 obs. of  3 variables:
 $ barcodes: chr  "AAACAGCCAAACTGCC-5" "AAACAGCCAACAGGAT-1" "AAACAGCCAAGATTCT-1" "AAACAGCCAAGGTATA-1" ...
 $ libcodes: Factor w/ 8 levels "1","2","3","4",..: 5 1 1 1 5 7
 $ samples : chr  "PN0252_0005" "PN0252_0001" "PN0252_0001" "PN0252_0001" ...
>

> pbmc
An object of class Seurat
36601 features across 26473 samples within 1 assay
Active assay: RNA (36601 features, 0 variable features)
>

> dim(metadata)
[1] 26473     3
>

Thanks @angel much appreciated

If you have single-dimension per-cell metadata, and it's arranged identically to the cell order in the Seurat object, I find it easier to use the double bracket notation to add metadata to a Seurat object. For example:
metadata$barcodes -> pbmc[["barcodes"]]
metadata$libcodes -> pbmc[["libcodes"]]
metadata$samples -> pbmc[["samples"]]

Answer from @user438383, converted from comment:
Did you try the AddMetaData function? This post may also be useful.
cluster_letters <- LETTERS[Idents(object = pbmc_small)]
names(cluster_letters) <- colnames(x = pbmc_small)
pbmc_small <- AddMetaData(
  object = pbmc_small,
  metadata = cluster_letters,
  col.name = 'letter.idents'
)


text:
Getting conventional gene symbol for Seurat
I have a Seurat object made by human single cells
I am mapping some genes on that but no sign of expression
When I GOOGLE for those genes I see the genes have different names
How I know Seurat uses which gene name please
This is my list
CD45
MHC II
CD11b
Ly6C
Ly6G
F4/80
CD11c
CD38
Arg1
SiglecF
CD206
CD62L
CD103
iNOS
PD-L1
TNFa
CD64
TCRgd
Foxp3
RORgt
CD8α
Tbet
CD25
IFN-γ
CD44
CD86
CD80
PD-1
B220
NK1.1
CD19
CD4
TCR β

From this list only 6 genes being mapped on the Seurat data

You have already asked this question and it has been answered - with the exact same response. Why are you opening multiple threads instead of exploring more on your own?

You might want to use Biomart to figure out the official names for those genes based on their aliases.

Can I trust NCBI? For instance this link which says TBX21 is an official gene name? https://www.ncbi.nlm.nih.gov/gene/30009

Yes, you can trust NCBI. Like I said in an earlier comment, HGNC provides the "official" gene symbols, and NCBI/EnsEMBL cross-referencing is reliable.


text:
Adding legend to Botton and top annotation in complexheatmap
I have a heat map but top and bottom annotation don't have legend
Top_ha = HeatmapAnnotation(Group = anno_simple(top$Group,
                                               col =structure(hue_pal()(length(unique(top$Group))),
                                                              names = unique(top$Group)),
                                           border = TRUE,
                                           height = unit(5, "mm")),
                       "Total Mutation Burden"=anno_barplot(clin$TMB, height = unit(15, "mm"), width=unit(30, "mm"),gp = gpar(fill = structure(hue_pal()(length(unique(clin$TMB))),
                                                                                names = unique(clin$TMB)))))

Bottom_ha = HeatmapAnnotation(Response = anno_simple(clin$Response,
                                                 border = T,
                                                 col = structure(c("purple4", "grey"), names = c("Responders", "NonResponders")),
                                                 height = unit(5, "mm"),show_),
                          "TP53 Status" = anno_simple(clin$TP53_status,
                                                          border = TRUE,
                                                          col = structure(c("firebrick", "grey"), names = c("Mutant","Wild type")),
                                                          height = unit(5, "mm")))

Right_Anno <- HeatmapAnnotation(Signature = anno_text(right$Signature), which = "row")

Heatmap(mat,
        bottom_annotation = Bottom_ha,
        top_annotation = Top_ha,
        name = "Exposure",
        right_annotation = Right_Anno,
        show_row_dend = FALSE,
        show_parent_dend_line = FALSE,
        column_split = c(top$Group),
        column_title = c("Mutagenic", "DDR Impaired", "C>A/T Dominant"),
        col = colorRamp2(c( 0, 0.5), c("white", "red")),
        show_column_names = FALSE)

As you are seeing Response, group and total mutation burden don't have any legend and annotation
Please can you help me to put legend for annotation?
I have tried annotation_lable but says Error: either provide levelsor provide namedcolors.
Also show_annotation_name = TRUE,show_legend = TRUE did not work

Answer from @ram-rs, converted from comment:
For the Group, try using it without the anno_simple. Just Top_ha = HeatmapAnnotation(Group=top$Group, "Total Mutation Burden" = ...), nothing else for Group. See if that generates a legend at all. We can customize later.hing else for Group. See if that generates a legend at all. We can customize later.


text:
Interpreting this plot from GSEA
I have RNA-seq for two groups of patients: Responders to chemotherapy (n=9) versus non-responders to chemotherapy (n=24)
I have done DESeq2 and obtain fold change of responders versus non-responders in which positive fold change means down-regulation in responders (if I am not wrong)
> head(df)
    gene       FC
1 286499 7.116071
2 260436 6.691189
3   2267 5.751322
4   5968 5.522079
5   4103 5.328008
6 442709 5.036087
>

I have done GSEA in pre-ranked way
got this plot for KRAS pathway

I interpret this as; KRAS pathway in being down regulated in responders
Am I right?
Please correct me
When I separated genes for up and down regulated
For up regulated genes I did not get any significant GSEA but for down regulated genes in responders I got this

SO does this mean that KRAS pathway is being up regulated in responders because I have used down regulated genes in responders?

A gene set enrichment analysis (GSEA) tests for enrichment of a gene set within a ranked list of genes.
The primary outcome of the analysis is enrichment or no enrichment. Gene sets frequently include correlated genes, and that correlation can be positive (enhancers or co-expressors) or negative (e.g repressors). It's reasonable to include negatively-correlated genes within gene sets because some pathways involve both enhancement of a particular biochemical pathway, and suppression of a pathway that carries out opposite things.
I don't know about the KRAS pathway you're comparing against here, but a statement I would make based on interpreting that GSEA graph would be something like, "The KRAS.600_UP.V1_DN gene set is enriched when comparing the expression of responders vs non-responders."
I'm not quite sure how you separated up and down-regulated genes (code would be helpful with that). Your down regulated gene plot looks like it's ranked in the opposite direction to the combined plot, which is not what I expect. I think that determining an enrichment score separately for upregulated and downregulated genes (within your ranked list of genes) is a good idea, but it shouldn't change the direction of the enrichment curve within each category. It will, however, alter the form of the curve, by stretching or compressing values vertically and altering the start/end offset.

In Deseq2, it depends on contrast eg. if your contrast <- c("Condition", firstC, SecondC) is then -ve is downregulated in FirstC and same for +ve. It is always better to take look at count data for the confirmation (responders versus non-responders samples).
[for more detail you can refer here][1]
[1]https://github.com/amarinderthind/RNA-seq-tutorial-for-gene-differential-expression-analysis/blob/master/RNA-Seq-DGE.R


text:
A GTF annotation for GRCH375d version
I have been given a bunch of .bam files, that come from the alignment on GRCh37d5 of human genome by Star tool. I want to extract raw read counts from these bam files by FEATURECOUNTS tool so I need a GTF. I googled a lot but I did not find GTF for GRCh37d5. I used another GTFs and I obtained weird features. Please one of you help me in finding this specific GTF. I converted a bam to fastq and I used Salmon to quantify gene expression by a transcripome.fasta but rate of mapping was 0.06

It's helpful to know that GRCh37 refers to the Genome Reference Consortium assembly of the human genome, version 37 (should be equivalent to hg19 from UCSC). I've never seen "5d" though. Often times GRC will introduce "patch" updates with minor changes that don't change the primary coordinate system—that's what p13 refers to in kayvee's answer below.

@DanielStandage Might be of interest to you: https://www.biostars.org/p/328824/#330791

The most recent annotation of GRCh37 assembly by NCBI RefSeq is here: https://ftp.ncbi.nlm.nih.gov/genomes/Homo_sapiens/GRCh37.p13_interim_annotation/
You can download the annotation in GFF3 format from that path and use it with featurecounts tool by first converting it to the SAF format using the commands shown below:
$ zgrep -v '^#' annotation.gff.gz \
    | awk 'BEGIN{FS="\t";OFS="\t"}($3=="exon"){print $1,$4,$5,$7,$9}' \
    | sed -r 's/ID.*GeneID:([0-9]*).*/\1/g' \
    | awk 'BEGIN{FS="\t";OFS="\t"}{print $5,$1,$2,$3,$4}'
    > annotation.saf

You should confirm that the chromosome sequence identifiers in the BAM file and the SAF file match. NCBI annotation has identifiers in the accession.version format whereas UCSC uses chr# format and Ensembl uses #. For this, you can use the -A parameter from featureCounts. Assuming that your BAM file uses UCSC-style identifiers, you can generate a 'chromosome name alias file' that can be used with featureCounts as follows:
## get NCBI assembly report for GRCh37
wget -nv ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.24_GRCh37.p12/GCF_000001405.24_GRCh37.p12_assembly_report.txt

## generate chr_alias.csv file
grep -v '^#' GCF_000001405.24_GRCh37.p12_assembly_report.txt \
    | dos2unix \
    | awk 'BEGIN{FS="\t";OFS=","}($10!="na"){print $10,$7}' \
    > chr_alias.csv

I didn't know about "SAF" format. I use gff with FeatureCounts. For the record, according to https://www.rdocumentation.org/packages/Rsubread/versions/1.22.2/topics/featureCounts: "The SAF annotation format has five required columns, including GeneID, Chr, Start, End and Strand. These columns can be in any order".


text:
Getting this information for GISTIC2
Long time I am struggling how to provide GISTIC2 required input like Num_Probes and Segment_Mean
I have called somatic copy number by scatngs and I have copy number output file (attached);
https://www.dropbox.com/s/r6nef2v1j2di6q6/my.copynumber.caveman.csv?dl=0
This is few lines of scatngs copy number output
1,1,13116,665401,2,1,3,1
2,1,668374,96324627,2,1,2,1
3,1,96326095,96419709,2,1,1,0

In scatngs reference paper says in this file we have
Segment number

Chromosome

Start position (origin-1)

End position (origin-1)

Major copy number—normal

Minor copy number—normal

Major copy number—tumor

Minor copy number—tumor

scatngs also returns another output like below
https://www.dropbox.com/s/7i28vooixdp06c2/WTSI-OESO_121_1pre.copynumber.txt?dl=0
Chromosome  Position    Log R   segmented LogR  BAF segmented BAF   Copy number Minor allele    Raw copy number
rs62635286  1   13116   -0.0290741822638953 0.245246748556562   0.25    NA  3   1   2.35806652809599
rs75454623  1   14930   0.402869048014209   0.245246748556562   0.3095  NA  3   1   4.14031506947036

I am wondering how I can get Num_Probes and Segment_Mean From this files?
People say
Segment_Mean = log2(tumour copy number / 2)

So, does that mean I should add up Major copy number—tumor and
Minor copy number—tumor to get tumour copy number? for getting Num_Probes what should I do?
EDITED

people says
scatngs uses read counts at SNPs to generate calls and so the
  locations of the SNPs is what I need to
  generate the number of markers in each segment. (eg each SNP = 1
  marker).
Do they mean, for instance if the position of SNP  rs62635286 is
  13116, so the num_prob is 13116?

Any help please?
Thanks a lot in advance

Cross-posted on GitHub for some reason: https://github.com/cancerit/ascatNgs/issues/84

I mentioned a cross-post. That does not stop anyone from reading the question or helping you.

Mentioning a cross post is hardly negative feedback or reason to ignore a post that belongs on bioinfo SE.

Did you try reading the GISTIC2 documentation about how to get the data in that format? Or perhaps they point to a known file extension (or a weird one from only a single program) I tried searching on a couple of website about these required Num_probes but couldn't find anything. Are you using this software or other with the same name?

I read GISTIC2 documentation, there is nothing about how to get required format. Exactly I am trying the same software you are referring to. People using GISTIC in broad institute produce their CN by this method that returns format for GISTIC https://software.broadinstitute.org/gatk/documentation/article?id=9143#1

Depends on the software to generate the segmentation file;
For example I have used scatngs for segmentation that for sample returns
Segment number
Chromosome
Start position (origin-1)
End position (origin-1)
Major copy number—normal
Minor copy number—normal
Major copy number—tumor
Minor copy number—tumor

By  log2(Major copy number—tumor +  Minor copy number—tumor /ploidy of each sample) you would get segment mean
For number of probes
The number of probes depends on where the data came from.
1) From SNP arrays:  Number of array probes within each segment.
2) From exon capture data:   Probes are the ends of exons and so the number of probes are equal to the number of exons  X 2.
3) From whole genome sequencing:  Probes are the boundaries of bins used to generate the copy number data and so the number of probes are equal to the number of bins +

For example in my case for scatngs
Each SNP is a marker, and so you have to count how many of these SNPs are on each segment

I have asked this question in
Extracting the number SNP in each range


text:
Extracting the number SNP in each range
I have called copy number and I have 2 files (I have shared the link of my files ); One contains some ranges
> head(cndata[,c(2,3,4)])
     Chromosome   Start      End
4470          1   51479   817980
4471          1  818499  1136753
4472          1 1138735  2558308
4473          1 2558740  5724264
4474          1 5724940  5749083
4475          1 5749226 12529544
>

https://www.dropbox.com/s/wei3i787c9nrz8q/WTSI-OESO_121_1pre.copynumber.caveman.txt.csv?dl=0
I have another file like below
              Chromosome         Position
rs62635286       1                 51479
rs75454623       1                 114930
rs806731         1                 30923

https://www.dropbox.com/s/tiqjhkrpfr9d9gu/WTSI-OESO_121_1pre.copynumber.txt?dl=0
In first  column I have some SNPs, how I can count the number of SNP in each range? For example based on the second table how many SNP are in range of 51479 to  817980?

Repost of https://www.biostars.org/p/376315/

We need to "merge with overlap/range" then "group by count", searching these terms should help you solve the problem.

Damian and Wouter have given you perfect solutions on BioStars. There is nothing more to add.

@FereshTeh I think people are getting a bit upset if you crosspost in different sites and it's not like you are a new user who just does not know. You are here for a while and you know how this site works. You must have noticed that crossposting is not very popular and if you do, at least you can explicitly mention it in the question.

Even better, if you get answers in one of the two places you ask, at least take the time to post them in the other. So if you get an answer here, also post it on your question at biostars and if you get an answer there, also post it here. That way, both communities benefit.

@FereshTeh I found a good article for you. Try always explain the problem (you do that moreless already), what have you tried to resolve it (this one you need to work on) and how to reproduce the problem you find (this one you do, but it would be nicer if you shared samples, not dropbox links).

There is like a thousand different ways how to achieve this. You could use a specialised software for this (like bedtools) or calculate it simply in R.
R solution: You can make a function that calculates the number of SNPs in a range and than you can apply it on all the ranges in the table with genomic ranges.
snp_table <- read.table('WTSI-OESO_121_1pre.copynumber.txt')
genomic_ranges <- read.table('WTSI-OESO_121_1pre.copynumber.caveman.txt.csv', sep = ',', col.names = c('Chromosome', 'Start', 'End'), stringsAsFactors = F)

how_many_in_range <- function(coords){
    coords = as.vector(coords)
    sum(snp_table$Chromosome == coords[1] & snp_table$Position > as.numeric(coords[2]) & snp_table$Position < as.numeric(coords[3]))
}

genomic_ranges$number_of_snps <- apply(genomic_ranges, 1, how_many_in_range)

Now the table genomic_ranges has a variable that sais how many SNPs are in the range.

You know, the problem others have with your question is not about your coding skills, but more about the crossposting that really annoys the heck out of everybody, also your questions do not clearly show much research effort. It might be that you spent days trying to figure it out, but from your question, we can not tell, because you did not write what you have tried. Finally, I don't think self-hatred is useful to improve your coding skills, reading manuals and trying it out is way more efficient.


text:
Showing name of genes in boxplot
I have a long format data frame like this
> head(df)
     gene

1:  CCL28   2.030674       2.023333 1.921508    1.952007
2:  CXCL2   2.204847       2.275296 2.175790    2.241554
3:  CXCR4   2.121186       2.147181 1.950065    1.982639
4: IFNGR1   2.311390       2.324111 2.220727    2.255556

And by this code
ggplot(df, aes(variable,value)) +
  stat_boxplot(geom="errorbar", width=.5)+
  geom_boxplot(aes(fill=variable))+
  theme_bw()+
  theme(axis.title.x=element_blank(), axis.title.y=element_blank())+
  stat_summary(fun.y=median, colour="red", geom="line", aes(group = 1))+
  geom_jitter(position = position_jitter(0.2))

I have this box
I have added this line to my code where I got
geom_text(data = df_melt, aes(x = variable, y = value, label = gene), size = 5)

[![enter image description here][2]][2]
I want to show the name of genes in box plot but I don't how, cam you help?
Thanks

I think you need to edit your question because the colnames of head(df) and the df you use in ggplot2 are not the same (no mention of the column variable and value. So, I believed that you reshape your dataframe. Please indicate all step and provide a reproducible example

In addition, to have several errors in your code for plotting (for example, you are using df for plotting and then df_melt for geom_text but you didnot mention what is it). Based on your graph, it seems that you have issues with the melting part.
However, these are quite secondary issues that are probably be fixed with the code below. with your example, there is one trick because you are using jitter to plot your points (that is not addressed by the link that @ATpoint provided in comments). This will randomly position the point at the right y value but with different x position in order to not overlap points.
So, when you will add your text label, the labeling won't fit to the position of points. But you can make the jitter reproducible by adding a seed (see this github discussion: https://github.com/tidyverse/ggplot2/pull/1996). However, geom_jitter function does not have the argument seed, so you can use geom_point function and the argument position_jitter instead.
So, using the few genes that you provided as an example. I have to reshape it (a step you probably did but you did not provide in your question), here I'm using pivot_longer function from tidyr package:
library(dplyr)
library(tidyr)
DF <- df %>% pivot_longer(., -gene, names_to = "variable", values_to = "value")

# A tibble: 16 x 3
   gene   variable       value
   <chr>  <chr>          <dbl>
 1 CCL28  Res_immune      2.03
 2 CCL28  Non_res_immune  2.02
 3 CCL28  Res_OBP         1.92
 4 CCL28  Non_res_OBP     1.95
 5 CXCL2  Res_immune      2.20
 6 CXCL2  Non_res_immune  2.28
 7 CXCL2  Res_OBP         2.18
 8 CXCL2  Non_res_OBP     2.24
 9 CXCR4  Res_immune      2.12
10 CXCR4  Non_res_immune  2.15
11 CXCR4  Res_OBP         1.95
12 CXCR4  Non_res_OBP     1.98
13 IFNGR1 Res_immune      2.31
14 IFNGR1 Non_res_immune  2.32
15 IFNGR1 Res_OBP         2.22
16 IFNGR1 Non_res_OBP     2.26

And then, for the plotting part, you can use geom_text_repel function of the ggrepel package to add labeling without overlapping your points:
library(ggplot2)
library(ggrepel)
ggplot(DF, aes(x = variable, y = value, fill = variable, label = gene))+
  geom_boxplot(alpha = 0.2)+
  geom_point(position = position_jitter(width = 0.2, seed = 2), show.legend = FALSE)+
  geom_text_repel(position = position_jitter(0.2, seed = 2), show.legend = FALSE)

Warning
If you have a long list of genes, it will be pointless to show all of them on the boxplot. It will makes things really illisible. If you want to show some particular genes of interest, make a subset of the original dataframe. A lot of posts on SO and here have addressed this subsetting.
Reproducible data example
structure(list(gene = c("CCL28", "CXCL2", "CXCR4", "IFNGR1"),
    Res_immune = c(2.030674, 2.204847, 2.121186, 2.31139), Non_res_immune = c(2.023333,
    2.275296, 2.147181, 2.324111), Res_OBP = c(1.921508, 2.17579,
    1.950065, 2.220727), Non_res_OBP = c(1.952007, 2.241554,
    1.982639, 2.255556)), row.names = c(NA, -4L), class = c("data.table",
"data.frame"), .internal.selfref = <pointer: 0x563d91b65350>)

Did OP add comments and then delete them? This seems to be a one-way conversation and OP either deleted their comments or edited their posts, adding their comments in-post.

You're right. OP deleted his comments. It was mainly asking for explanation about error message that I presumed have been solved. I will remove my comments and they are not anymore relevant

That behavior is both entitled and ungrateful. Wow.


text:
Extracting some part of a list
I have two normalized gene expression values (log2 of cpm)
> head(biomarker[,1:5])
               A2        A3        A4        A6        A7
A2M    12.4071618 12.631601 12.889748 11.842511 11.574134
ABCB11  0.5151957  6.176788  5.414905  7.134915  1.247590
ABCC2   6.7244303  6.487794 10.178615  6.384132  6.089462
ABCC6   6.3977122  5.823627  7.370761  7.397091  6.071587
ABCF1   9.1609847  9.893258 10.116638 10.225520 10.838486
ABCG2   7.6841874  7.738293  6.833152  6.888041  6.030032
>

And
> head(immune[,1:5])
              A2        A3        A4        A6        A7
A2M    10.904748 11.388404  9.910614 11.513439 12.963609
ABCB11  5.011380  6.359443  7.145992  8.451947  7.722605
ABCC2   5.040461  6.477014  3.873996  6.409777  9.133971
ABCC6   7.798441  7.601848  9.948072 11.628533 12.701460
ABCF1   8.553597  8.615425 11.145903 10.289098 11.444140
ABCG2   6.224294  5.629375  8.293416  7.979859  8.603793
>

> dim(biomarker)
[1] 719  56
> dim(immune)
[1] 719  56
>

I have matched samples in these experiments and same gene too. I could roughly say in both experiments everything has been constant. However, I want to do Student's t-Test to find genes that they had inconsistent gene expression patterns between the two data sets.
By this function I am getting a list but I don't know how to extract the names of genes with p-value > 0.05 or whatever cut-off
> f <- function(x,y){
+     test <- t.test(x,y, paired=TRUE)
+     out <- data.frame(stat = test$statistic,
+                       df   = test$parameter,
+                       pval = test$p.value,
+                       conl = test$conf.int[1],
+                       conh = test$conf.int[2]
+     )
+     return(out)
+ }
> t=sapply(seq(ncol(t(biomarker))), function(x) f(t(biomarker)[,x], t(immune)[,x]))

Any help please?
  str(t)

    [list output truncated]
     - attr(*, "dim")= int [1:2] 5 719
     - attr(*, "dimnames")=List of 2
      ..$ : chr [1:5] "stat" "df" "pval" "conl" ...
  ..$ : NULL
    >

> sapply(t, "[[", "pval")
Error in FUN(X[[i]], ...) : subscript out of bounds
>

> pvalueGenes <- sapply(colnames(biomarker),function(i) t.test(biomarker[i, ],immune[i, ], paired = TRUE)$p.value)
Error in if (stderr < 10 * .Machine$double.eps * abs(mx)) stop("data are essentially constant") :
  missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In mean.default(x) : argument is not numeric or logical: returning NA
2: In if (stderr < 10 * .Machine$double.eps * abs(mx)) stop("data are essentially constant") :

 Show Traceback

 Rerun with Debug
 Error in if (stderr < 10 * .Machine$double.eps * abs(mx)) stop("data are essentially constant") :
  missing value where TRUE/FALSE needed

Did you try iterate by colnames instead? Have you searched the internet about how to get names of the variables you iterate over them? BTW you might be interested in the broom package, function tidy.

Sorry I don't know how to do that

A simple google search would give you answers. You could use tibble::rownames_to_column(biomarker, col_name). You could even do biomarker$col_name <- colnames(biomarker)

Sorry what is col_name please?

@FereshTeh you don't need to apologise for not knowing something. "col_name" is an abbreviation for column name. Use the help pages of those functions Ram mentioned to understand what do you need and what do they do

@FereshTeh Did you read and understood the help page of the rownames_to_column function?

Actually, no; You know, likely I have done t-test to find differentially expressed genes. I have 719 genes and two sets of samples while the function returns a list with 56 * 4 dimension that means t-test has been done between samples while I need to extract genes with p-value > 0.05 among two sets of datasets, that is why I got confused and I don't know rownames_to_colnames would solve which of these

@FereshTeh Ok, that is something you should [edit] into the question.

If you're going to downvote, add a reason why. Don't be passive-aggressive. It's bad for the community. Thanks!

@AlexReynolds It is not required to disclose, but I think that the downvote button already explain this: "does not show research effort, it is unclear or not useful". Asking for an explanation from downvoters is also passive-aggressive...

@FereshTeh I used col_name as a name for the new column that will be generated by the function rownames_to_column. I should have used 'col_name' (with quotes), sorry about that!

@AlexReynolds I kinda agree with llrs. In addition to not showing research effort, OP has a habit of editing her posts extensively and changing the entire premise once she receives a response, and such behavior needs to be discouraged.

Starting from the definition of t:
...
> t <- sapply(seq(ncol(t(biomarker))), function(x) f(t(biomarker)[,x], t(immune)[,x]))

Set the column names to the gene names from one or the other of your conditions:
> colnames(t) <- rownames(immune)

Set the cutoff:
> cutoff <- 0.05

Filter for rows which meet the p-value cutoff criterium, and get their column names (gene names):
> colnames(t)[which(unlist(lapply(t["pval",], function(x) x < cutoff)))]
[1] "ABCC6"
...

Only one gene name is reported here, using your demo datasets. Adjust your cutoff and criterium logic, as needed. For example, I filter here for gene names where the p-value is less than the cutoff; however, from your question, you may be asking for the inverse, for whatever reason, so you could change < to >=.

I might suggest changing the variable name t to something else, as t() takes the transpose of the argument, if of matrix or data frame type. Using non-ambiguous variable names will make it easier for you to make changes or fix problems down the road.

From your example we are comparing samples not genes, I think we need to run the test row wise, something like this:
pvalueGenes <- sapply(colnames(biomarker),
                      function(i) t.test(biomarker[i, ],
                                         immune[i, ], paired = TRUE)$p.value)

names(pvalueGenes)[ pvalueGenes < 0.05 ]

Thank you, I added output of what you suggested in my post

Sorry, I have 2 matched matrices with matched samples and genes I want to find genes with p-value < 0.05 between 2 matrices but I could not do that

@FereshTeh forget your function f, there is no need, just run the commands in my answer. And check what you see in pvalueGenes object.

Thank you but your code itself giving error please have a look at error I added to my post

@FereshTeh The reason for error is statistics, read about it more here. Suggestion, when you get an error try to use the search, 99% of the time the error has been searched before and you will find the solution/explanation.

@zx8754, although names(x)[ x < 0.05 ] works, would it not be better to use names(x[ x < 0.05 ]), so each layer is a logical unit?

@RamRS Not sure what you mean by "logical unit". Yes both approaches work, but your version is about 1.14 times slower, tested on 100K vector.

Ah, I see. Never mind, then. By logical unit, I meant that x, x[x<0.05] and names(obj) are all self-contained, whereas names(x)[boolean_vector] is written in a way that the subsetting depends on understanding the statement as a whole.


text:
What does that mean if there is no bifurcations in my time series single cell data?
Sorry,
I have analysed my time series single cell data with Scuba algorithm but that says that there is no bifurcations in my data and outputted a single chain of 4 states with no branching like below figure

Sorry, I googled but I am not quite sure what this means
Any suggestion please?
Thanks a lot

The most obvious interpretation is that the four cell lines are related in a linear pathway rather than a branching one.


text:
Subset cells based on a changefold level in gene expression
is there a way to subset a group of cells based on the foldchange of a gene expression? I have the KRAS object KRAS An object of class Seurat 53805 features across 6826 samples within 1 assay Active assay: RNA (53805 features). I want to subset on the expression of the Olfm4 gene, but with a statistical threshold, for example a logfold2 change of 1. In the meantime I use subset(KRAS,subset='Olfm4'>1).
Thanks

Please [edit] your question and give us a few more details. What data do you have? In what format? You tagged with seurat, so I guess you are using that tool. If so, please tell us! Ideally, show us a minimal, reproducible example of your code and data so we can test our answers and make sure we get what you need.

Could you post how do you have your data which column has this fold change and what what do you expect to achieve? Also it would be helpful to think about the fold change, what are you comparing?

Hi @IroPierides, Seurat is quite an easy question for this board. If you detail the experiment and/or supply data, I'm sure they'll help. Note very similar questions have been answered in past posts.


text:
No output when trying to obtain protein sequence from PDB file
I am trying to obtain the sequence from PDB file, but no output is given when I run the following code :
from Bio import SeqIO

with open("1a9l.pdb", "r") as handle:
    for record in SeqIO.parse(handle, "pdb-atom"):
        print(record.seq)

If you only want protein sequences for an PDB code, NCBI has already a dataset (in genbank format) for that. Therein the accession IDs are the code underscore chain letter...

I was trying to do it through the script,as it is part of pipeline..

If you are using a conda version of python you can use Pymol module for this, which is actually more robust in some aspects than bio.PDB, although errors may cause segmentation faults. The following:
import pymol2
## paralellisable version.
with pymol2.PyMOL() as pymol:
    pymol.cmd.fetch('1a9l')
    print(pymol.cmd.get_fastastr('polymer.nucleic'))

Gives me:
>1a9l_A
GGGUGACUCCAGAGGUCGAGAGACCGGAGAUAUCACCC

Alternatively, you can use PDBe's API.
import requests

response = requests.get('https://www.ebi.ac.uk/pdbe/api/pdb/entry/molecules/1a9l').json()['1a9l']

for entity in response:
    if entity['molecule_type'] == 'polyribonucleotide':
        print(entity['sequence'])

Do note that missing 'residues' (residues, small molecules and nucleotides) will cause issues. PDB files do not have missing 'residues', while mmCIF files have a dictionary that does. The API will give you the latter.

Hi, I tried conda install -c schrodinger pymol to install pymol module. On import, it shows me an error : ModuleNotFoundError: No module named 'pymol2'
Is this the way to download it for use with conda?

conda install -c schrodinger pymol will install both. Does import pymol work? That is, is the version of python you are running the conda one? If so, single threaded pymol in the example runs as import pymol without the context manager line. Doing conda search -c schrodinger pymol for me (Ubuntu) gives me choices including pymol 2.1.1 for  py27_2 and py36_2. Do you get the same versions?

yes,I get same.No import pymol is also not working,yes I am running on same version of python as conda.

Okay, just to clarify: (a) there is no sudo'ing/nohup/systemctl service going on, (b) your which python (or whichever you are calling) and your which conda point to the same folder and (c) when you type conda list pymol it gives it as installed, correct?

I should have phrased it not as which python, but within your run sys.executable gives that one.

yes,it shows:
Name                    Version                   Build  Channel
pymol                     2.3.4            py37h95b93ae_0    schrodinger

Installed with conda pymol is no different than numpy or pandas, so this is odd if the python is indeed correct one (cf. sys.executable ) and I have never had issues with it on dozens of systems (MacOS, Debian, RedHat) and architectures (x64, arm64) w/o issues. Uninstalling and reinstalling it?

Does pip freeze | grep 'pymol' - give the same result? Does print([k for k in sys.modules.keys() if 'pymol' in k.lower()]) gives you something? If so, pymol2 = importlib.import_module('pymol2') just might work, but is nearly semantically equivalent to import pymol2 so I highly doubt it.

I will try to troubleshoot this issue this week.If import issue still persists,will ask you.

That PDB entry is an RNA molecule, so there is no protein sequence to extract.
If it extracted the RNA sequence by default here you can see there could be an ambiguity as to whether it was an RNA sequence or a protein sequence, as the alphabets share characters.

Got it..So,I can't extract RNA sequence through this option.

Not using SeqIO from Biopython, no. The "pdb-seqres" option won't work here either. You could download the flat file of PDB sequences and read them into a Python dictionary: https://www.rcsb.org/pages/download_features#FASTA. Be aware though that sometimes atoms are missing in the ATOM records, which may have been why you wanted to read directly from the atoms in the first place.


text:
How are different forcefields developed/evaluated?
The Rosetta software can use a variety of different energy functions:

Ref2015
GenPot
Beta-Nov16

How are these energy functions developed and validated? Is it common to do a single protein design with several of these energy functions?

To misquote Bismarck, "If you like forcefields and sausages, you should never watch either one being made".

Some force-fields rely solely on theoretical physical principles, while others are tweaked specifically for the goal of accurately depicting certain structures like sheets and helices or movement/diffusion aspects at different theoretical temperatures of a set of well-studied proteins. (Top-Down vs. Bottom-Up)
I haven't done MD myself but I am in an MD group and people seem to normally stick to a single force-field for a given MD project but it happens that they do very short simulation runs with force-fields just as a test.
They definitely don't plan on doing the full simulation on different force-fields unless the comparison of those is their main goal.

The Rosetta ref2015 is well described in two papers, which I wish I had read much earlier, see this JChIM paper.
It's a statistical forcefield, where the terms are weighted to best match the empirical total energy in kcal/mol and the paper does a swell job at explaining what they do.
For several Rosetta experiments, you switch between full-atom and centroid (=coarse grain). These have different weights, as you can see from the files in the database folder or via pyrosetta, but majorly they have different "params" (topologies) for the residues as the centroid lacks many atoms as it's a simplification. In some instances, you may switch between dihedral (normal) to cartesian (explained in the FastRelax paper), wherein the former used internal coordinates (e.g. atom A is 1.5 Å away from B and A,B,C is a 107ª angle etc —cf. a params file ICOOR block), while the later is cartesian coordinates ( A is at [x, y, z] —cf. PDB file).
The GenPot uses its own AtomTypes, which means it uses its own params for each residue —like centroid/fa. The topology of a residue (ResidueType) is based on its conformation (ICOOR), bonding and atom definition, which is atom name, gasteiger-marsili charge and AtomType (an amalgamation of element, sp-hybridization, size, H-bond donor/acceptor etc.) These AtomTypes are better for ligands as many groups behave different from peptide atoms —say, a nitrile nitrogen is definitely different from a deprotonated aromatic nitrogen.
Some application, such as hydration, ligand docking etc. have specifically tweaked scorefunctions for use by their movers. FastRelax is a nice example as the cycle weighting changes are spat out by the tracer at vertiginous speed, where the protocol is:
repeat 5
coord_cst_weight 1.0
scale:fa_rep 0.040
repack
scale:fa_rep 0.051
min 0.01
coord_cst_weight 0.5
scale:fa_rep 0.265
repack
scale:fa_rep 0.280
min 0.01
coord_cst_weight 0.0
scale:fa_rep 0.559
repack
scale:fa_rep 0.581
min 0.01
coord_cst_weight 0.0
scale:fa_rep 1
repack
min 0.00001
accept_to_best
endrepeat

As you can see, the repulsion term of the Lenard-Jones is going from 0.04 up. These values were chosen not for thermodynamic reasons, but worked best for the application (energy minimisation of the structure).
However, this switching is all internal. Scores are always reported with the same scorefunction and comparing between different ones is not advised even for the minority of scorefunctions calibrated to be in kcal/mol.


text:
How to find a homolog that has a PDB structure available
I tried to find a homolog that has a PDB structure available so I can use this PDB file for comparative modeling. I ran a BLAST search but none of the search results seemed to have a known structure. Is there any way to find the results I want?

If you are using NCBI Blast and blastp against the pdb database gives you nothing, you can do a psiblast against the normal refseq database, do a few rounds, save the PSSM matrix and submit a new job with the PSSM matrix but against the PDB database.
You might need to limit the species of you get lots of too close homologues —say blasting E. coli protein for which you'd want to exclude Enterobacteriales.

You could use HMM-based comparison, for example with HHpred. This is generally the best way to find structural templates because structure is more conserved than sequence, and comparing HMMs will turn up hits of the same structure but minimal sequence overlap.
If HHpred turns up nothing you are in the region of de novo modelling, which is a different problem. However, it has advanced in recent years to the point where it may be useful to you. For example, you could try our method DMPfold.

To expand the first answer by @MatteoFurla: re-run the P-blast and under database option, which is set by default to non-redundant database (nr) select Protein databank (PDB), the third option from the bottom of the menu. All the Blast results will be from PDB and thus all results will have associated structures. Also note that the Genbank codes will be PDB codes.
If you have zero hits then consider PSI-Blast, how I suspect you are running P-Blast via the nr database.


text:
How to filter the PDB for single-domain proteins with full length 3D-structures solved
I am trying to reproduce a machine learning model that has been developed here. As one of the datasets, they use single-domain proteins with full length 3D-structures solved. Since I don't have a background in biology, I am not sure how to filter the PDB for these kind of proteins. Any help would be appreciated.

If you want to reproduce the model in that paper, you should email the corresponding author asking for the exact set of PDB IDs/chain IDs/residue ranges they used. There is too much flexibility in the criteria as given to get the exact set any other way.
If you want to assemble a similar set for your own work, the RCSB PDB advanced search is very powerful and would let you filter e.g. by max sequence length or for only monomers.
However the PDB doesn't deal with domains (regions that evolve and fold independently) directly. A PDB chain may contain part of a domain, one domain, or multiple domains.
If you want a dataset of protein domains you could use CATH, SCOP or ECOD.
As a side issue, please remember to remove overlap of homologous proteins between the training and test sets. People entering the field often neglect to do that. See for example this comment piece.

That is very helpful. Thank you for your answer.