Dear Friederike,

thanks for your reply. When I wanted to post my code here, I found a bad mistake in my code that change all my results, so I've edited the second figure. However, my previous question still remains: which normalization is correct? here is my code:

First strategy

query <- GDCquery(project = "TCGA-COAD" ,
                      data.category = "Transcriptome Profiling" ,
                      data.type = "Gene Expression Quantification",
                      workflow.type = "HTSeq - Counts" ,
                      sample.type = c("Primary Tumor", "Solid Tissue Normal"), 
                      experimental.strategy = "RNA-Seq")

    GDCdownload(query, method = "api")

    query.counts.colon <- GDCprepare(query)

    Colon.Matrix <- as.data.frame(SummarizedExperiment::assay(query.counts.colon))

    #sample grouping 
    dataSmTP <- TCGAquery_SampleTypes(getResults(query,cols="cases"),"TP")
    dataSmNT  <- TCGAquery_SampleTypes(getResults(query,cols="cases"),"NT") 

    #building an ordered matrix based on grouping 
    NormalGroup <- Colon.Matrix[,dataSmNT]
    TumorGroup <- Colon.Matrix[,dataSmTP]
    Sorted.Matrix <- cbind(NormalGroup,TumorGroup)

    #Normalization by DESeq2
    gr <- factor(c(rep("Normal", length(dataSmNT)) , rep("Tumor", length(dataSmTP))))
    colData <- data.frame(group=gr, type="paired-end")
    cds <- DESeqDataSetFromMatrix(Sorted.Matrix, colData, design= ~group)
    cds <- DESeq(cds)
    cnt <- log2(1+counts(cds, normalized=TRUE))

    dif <- data.frame(results(cds, c("group", "Normal" , "Tumor")))
    dif$padj <- p.adjust(dif$pvalue, method = "BH")
    dif <- dif[order(dif$padj),]


    #volcano plot
    ggplot(dif, aes(log2FoldChange, -log10(pvalue) , color=log2FoldChange))+geom_point()+theme_bw()

Second:

#Here is my code for matched samples
query <- GDCquery(project = "TCGA-COAD" ,
                  data.category = "Transcriptome Profiling" ,
                  data.type = "Gene Expression Quantification",
                  workflow.type = "HTSeq - Counts" ,
                  sample.type = c("Primary Tumor", "Solid Tissue Normal"), 
                  experimental.strategy = "RNA-Seq")

GDCdownload(query, method = "api")

query.counts.colon <- GDCprepare(query)

Colon.Matrix <- as.data.frame(SummarizedExperiment::assay(query.counts.colon ))

#sample grouping 
dataSmTP <- TCGAquery_SampleTypes(getResults(query,cols="cases"),"TP")
dataSmNT  <- TCGAquery_SampleTypes(getResults(query,cols="cases"),"NT") 


#comparing barcodes for normal and tumor samples
#separating the matched samples based on TCGA barcode
mached_tumor <- c()
for (i in 1:length(dataSmNT)){
  split_normal <- unlist(strsplit(dataSmNT [i],"-"))
  for(j in 1:length(dataSmTP)){
    split_tumor <- unlist(strsplit(dataSmTP [j],"-"))
    if(split_normal [2] == split_tumor[2] & split_normal [3] == split_tumor[3] &
       split_normal [6] == split_tumor[6] & split_normal [7] == split_tumor[7]){
      mached_tumor <- c(mached_tumor, dataSmTP[j])
    }
  }
}

#building an ordered matrix based on grouping 
NormalGroup <- Colon.Matrix[,dataSmNT]
TumorGroup <- Colon.Matrix[,mached_tumor]
Sorted.Matrix <- cbind(NormalGroup,TumorGroup)


#Normalization by DESeq2
gr <- factor(c(rep("Normal", length(dataSmNT)) , rep("Tumor", length(mached_tumor))))
colData <- data.frame(group=gr, type="paired-end")
cds <- DESeqDataSetFromMatrix(Sorted.Matrix, colData, design= ~group)
cds <- DESeq(cds)
cnt <- log2(1+counts(cds, normalized=TRUE))

dif <- data.frame(results(cds, c("group", "Normal" , "Tumor")))
dif$padj <- p.adjust(dif$pvalue, method = "BH")
dif <- dif[order(dif$padj),]

#volcano plot
ggplot(dif, aes(log2FoldChange, -log10(pvalue) , color=log2FoldChange))+geom_point()+theme_bw()

You are right, they shouldn't be the same, but indeed I don't know which one is correct.