Hi All! I need a bit of help getting my head around this problem.
We did TERF2 gene silencing (3 biological replicates, both control case ie 3 sets) in HT1080 fibrosarcoma cell line using shRNA that binds with full complementarity to its 7th exon. In WB we see 30-40% protein level downregulation (in two sets) and 20-30% in (3rd set). qPCR also shows 60% silencing in 2 sets and 40% silencing in one set.
However, when I performed Differential Expression analysis, I couldn't see TERF2 downregulated significantly. My log2FC was -0.34249 with a p-value of 0.075 and p-adj 0.721. When I compared raw read counts the difference between control &case in all three replicates was about 200-300 reads. What could have gone wrong? Out data is paired end RNAseq
We suspected that there might be a novel transcript that has not been reported so far and isn't degraded and might be contributing to increased read numbers. For that, we checked the Isoform percentage values that RSEM produces. We found that majority of the transcripts were coming from.
I think TERF2-201 is going down while TERF-2-206 is coming up, which is visible in two sets, but I am not sure.
In addition to the questions above, I would like to know if the majority of reads (~70%) are coming from a single transcript, is it true that if we count reads exon-wise (within exon boundaries), the number of reads (or distribution) will be similar across all exons.
rsem-prepare-reference --gtf /path to/gencode.v33.annotation.gtf --star --star-path ~/anaconda3/bin/ /path to/GRCh38.primary_assembly.genome.fa /home/parashar/scratch/name_of_output_file 2> rsem.stderr
rsem-calculate-expression --star-output-genome-bam --strandedness reverse -p 32 --calc-pme --calc-ci --keep-intermediate-files --append-names --sort-bam-by-coordinate --paired-end --estimate-rspd --star --star-path ~/anaconda3/bin/ ~/trimmomatic/paired/Contr1_S15_L004_R1_p.fastq ~/trimmomatic/paired/Contr1_S15_L004_R2_p.fastq ~/rnaseq/rsem ~/Contr1_S15_L004 2> ~/Contr1_S15_L004.stderr
The DESeq codes are here
Please Note the DE-Seq didn't account for batch correction.I am on it.
The HTSeq-count was used for raw feature estimation in union mode.
The RSEM was used in STAR mode.
The Star genomic alignment was used for downstream analysis
Alignment Stats from STAR Aligner
Obviously a batch effect as you say, try to remove it with removeBatchEffect on the vst counts and plot the PCA again, see if this removes it. If so then you can include it into the model as
~ batch + YourCondition.As indicated by ATpoint, there seems to be a batch effect. Also note that the 7th exon is not present in all transcripts and that several transcripts are annotated as being incomplete. So you actually don't know the full set of transcripts you're hitting with your shRNA. Add on top that the human reference genome is not the genome of the cell line you're using so there may be some differences. If you want to sort this out properly, you'll need to figure out the full sequence of all transcripts expressed in your cell line under your culture conditions. But is it worth the trouble? I would suggest to take a step back. Does it really matter in the context of the experiment? I mean what biological question are you trying to address and what do you care about in this context? If the gene function is mediated by proteins AND your antibody detects all relevant isoforms (some of the transcripts are not supposed to produce any protein and the function may only be carried out by some known isoforms) then you can assume you've removed ~30% of the function (assuming no complications such as post-translational modifications). That might be all you need to interpret your experiment.
We want to study which genes are differentially upregulated on TERF2 Silencing. But before that, it is essential to confirm TERF2 silencing happened or not. Moreover, when I run HTSeq, isn't it blind to isoform information and simply counts the PE primary alignments falling within the genic region. I run it in union mode.
What possible informatics approach do you suggest to probe into the cause? I am performing the batch correction, but do you think a mere difference of 200-300 raw read counts (between control-case) would fill up for this gene and it will be seen downregulated?
You see it on Western blot, so why bothering with RNA-seq? You could go for differential isoform expression with DEXseq though.
Yes I am on it as well.
Typically for protein-coding genes, silencing means reducing the amount of protein produced because we're interested in the function carried out by the proteins. You seem to have confirmed that some of the protein is gone. For me that's enough unless you care about the non-coding products of the genes. RNA seq is just taking you down a rabbit hole. As I wrote above, you have several unknowns. Take another genome annotation as reference and you may get a different result. If you insist in reconciling the Western blot with your RNA seq data, you need to know which transcripts are expressed in your cell line under your experimental conditions and verify that the treatment doesn't affect your target gene expression. If you really want to be sure of the Westen blot result, do a rescue experiment: mutate the gene to make it resistant to your shRNA and see that now you don't get a reduction in the protein.
Yes I have suggested the same experiment; to use Luciferase Assay and see if in the presence of shRNA the luciferase expression is repressed or not compared to wild type. However, I learned from them that they have done it multiple times and even have created stable cell lines where TERF2 is silenced.