|
转录组如果只看表达量真的是超级简单，真是超级简单，而且人家作者本来就测是SE50，这种破数据，也就是看表达量用的！ 首先作者分析结果是： 数据在GEO地址是：https://www.ncbi.nlm./geo/query/acc.cgi?acc=GSE50177 我们需要下载的RNA-seq的数据： https://www.ncbi.nlm.//sra/?term=SRP029245 https://trace.ncbi.nlm./Traces/study/?acc=SRP029245 ftp://ftp-trace.ncbi.nlm./sra/sra-instant/reads/ByStudy/sra/SRP/SRP029/SRP029245 下载地址很容易获取啦！ for ((i=677;i<=680;i++)) ;do wget ftp://ftp-trace.ncbi.nlm./sra/sra-instant/reads/ByStudy/sra/SRP/SRP029/SRP029245/SRR957$i/SRR957$i.sra;done ls *sra |while read id; do ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump –split-3 $id;done 因为我用fastqc看了看数据质量，发现没有什么问题，代码如下： ls *fastq |xargs ~/biosoft/fastqc/FastQC/fastqc -t 10 所以直接用hisat2软件把测序得到的fastq文件比对到hg19参考基因组上面 reference=/home/jianmingzeng/reference/index/hisat/hg19/genome ~/biosoft/HISAT/current/hisat2 -p 5 -x $reference -U SRR957677.fastq -S control_1.sam 2>control_1.log ~/biosoft/HISAT/current/hisat2 -p 5 -x $reference -U SRR957678.fastq -S control_2.sam 2>control_2.log ~/biosoft/HISAT/current/hisat2 -p 5 -x $reference -U SRR957679.fastq -S siSUZ12_1.sam 2>siSUZ12_1.log ~/biosoft/HISAT/current/hisat2 -p 5 -x $reference -U SRR957680.fastq -S siSUZ12_2.sam 2>siSUZ12_2.log 而且查看log日志可以发现，比对效果杠杠的： 93.10% overall alignment rate 然后把sam文件根据reads name来排序并且转换为bam文件节省空间 ls *sam |while read id;do (nohup samtools sort -n -@ 5 -o ${id%%.*}.Nsort.bam $id &);done 最后用htseq-counts工具来对每一个样本进行基因的表达量定量！ ls *.Nsort.bam |while read id;do (nohup samtools view $id | ~/.local/bin/htseq-count -f sam -s no -i gene_name - ~/reference/gtf/gencode/gencode.v25lift37.annotation.gtf 1>${id%%.*}.geneCounts 2>${id%%.*}.HTseq.log&);done 得到的文件如下： 这4个样本的基因的counts数据就可以用一系列的R包来做差异分析了，包括limma的voom，DEseq2，edgeR等等。这些包的用法都烂大街了，我就不赘述了。 做完差异分析，就可以跟作者的结果做对比，看看自己做的是不是对的。 其实你还会面临一些小问题，比如作者给的差异基因是用refseq ID表示的呀，如果注释成基因的symbol，还有name呢？ |
|
|
