comparison report_clonality/RScript.r @ 69:a59de79f6c0f draft

Uploaded
author davidvanzessen
date Mon, 08 Aug 2016 08:53:33 -0400
parents 69a6bba9ada9
children 7c9b9d97d596
comparison
equal deleted inserted replaced
68:85a1d832d197 69:a59de79f6c0f
570 } 570 }
571 571
572 res[is.na(res)] = 0 572 res[is.na(res)] = 0
573 infer.result = infer.clonality(as.matrix(res[,2:ncol(res)])) 573 infer.result = infer.clonality(as.matrix(res[,2:ncol(res)]))
574 574
575 print(infer.result)
576
575 write.table(data.table(infer.result[[12]]), file=paste("lymphclon_clonality_", sample_id, ".csv", sep=""), sep=",",quote=F,row.names=F,col.names=F) 577 write.table(data.table(infer.result[[12]]), file=paste("lymphclon_clonality_", sample_id, ".csv", sep=""), sep=",",quote=F,row.names=F,col.names=F)
576 578
577 res$type = rowSums(res[,2:ncol(res)]) 579 res$type = rowSums(res[,2:ncol(res)])
578 580
579 coincidence.table = data.frame(table(res$type)) 581 coincidence.table = data.frame(table(res$type))
773 write.table(newData, "junctionAnalysisUnProd_median.csv" , sep=",",quote=F,na="-",row.names=F,col.names=F) 775 write.table(newData, "junctionAnalysisUnProd_median.csv" , sep=",",quote=F,na="-",row.names=F,col.names=F)
774 } 776 }
775 777
776 PRODF = bak 778 PRODF = bak
777 779
780
781 # ---------------------- D reading frame ----------------------
782
783 D.REGION.reading.frame = PRODF$D.REGION.reading.frame
784
785 D.REGION.reading.frame[is.na(D.REGION.reading.frame)] = "No D"
786
787 D.REGION.reading.frame = data.frame(table(D.REGION.reading.frame))
788
789 write.table(D.REGION.reading.frame, "DReadingFrame.csv" , sep="\t",quote=F,row.names=F,col.names=T)
790
791 D.REGION.reading.frame = ggplot(D.REGION.reading.frame)
792 D.REGION.reading.frame = D.REGION.reading.frame + geom_bar(aes( x = D.REGION.reading.frame, y = Freq), stat='identity', position='dodge' ) + ggtitle("D reading frame") + xlab("Frequency") + ylab("Frame")
793
794 png("DReadingFrame.png")
795 D.REGION.reading.frame
796 dev.off()
797
798
799
800
778 # ---------------------- AA composition in CDR3 ---------------------- 801 # ---------------------- AA composition in CDR3 ----------------------
779 802
780 AACDR3 = PRODF[,c("Sample", "CDR3.Seq")] 803 AACDR3 = PRODF[,c("Sample", "CDR3.Seq")]
781 804
782 TotalPerSample = data.frame(data.table(AACDR3)[, list(total=sum(nchar(as.character(.SD$CDR3.Seq)))), by=Sample]) 805 TotalPerSample = data.frame(data.table(AACDR3)[, list(total=sum(nchar(as.character(.SD$CDR3.Seq)))), by=Sample])