链特异性文库
数据对应的论文
Identification and characterization of novel lncRNAs in Arabidopsis thaliana
看了下投稿时间
image.png
3天，也太吓人了吧！！！
期刊是 Biochemical and Biophysical Research Communications
BBRC 印象里好像看过介绍这个期刊的文章，特点就是审稿快，4区，3年平均影响因子2.577，版面费挺贵，2000多美元。
以上是题外话
下面开始数据分析流程
下载测序数据
~/.aspera/connect/bin/ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR547/009/SRR5478559/SRR5478559_1.fastq.gz ./
~/.aspera/connect/bin/ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR547/009/SRR5478559/SRR5478559_2.fastq.gz ./
NCBI下载拟南芥参考基因组
为了减小数据量，只用拟南芥的第一条染色体的前40000个碱基
准备gff3格式的注释文件
grep "CP002684.1" TAIR10.1_genomic.gff > chr1.gff
head -n 213 chr1.gff > chr1_pra.gff
可视化展示gff3文件
library(GenomicFeatures)
txdb<-makeTxDbFromGFF(file="chr1_pra.gff",format="gff3")
library(ggbio)
autoplot(txdb,
which=GRanges("CP002684.1", IRanges(100, 40000)),
names.expr = "gene_id",aes(fill=gene_id))+
theme_bw()
image.png
准备参考序列
samtools faidx TAIR10.1_genomic.fna
samtools faidx TAIR10.1_genomic.fna CP002684.1:1-40000 > chr1_pra.fna
参考基因组构建索引
hisat2-build Reference/chr1_pra.fna Reference/chr1
mkdir output
使用fastp对数据过滤
fastp -i raw_data/SRR5478559_1.fastq.gz -I raw_data/SRR5478559_2.fastq.gz -o raw_data/SRR5478559_clean_R1.fastq -O raw_data/SRR5478559_clean_R2.fastq
比对
hisat2 -p 12 -x Reference/chr1 -S output/SRR5478599.sam -q --no-unal --dta --rna-strandness RF -1 raw_data/SRR5478559_clean_R1.fastq -2 raw_data/SRR5478559_clean_R2.fastq
samtools view -S -b -o output/SRR5478599.bam output/SRR5478599.sam
samtools sort -O bam -o output/SRR5478599.sorted.bam output/SRR5478599.bam
R语言可视化bam文件
library(GenomicAlignments)
ga<-readGAlignments("SRR5478599.sorted.bam",use.names=T)
library(ggbio)
autoplot(ga,aes(fill=strand))+
theme_bw()
image.png
stringTie组装转录本
stringtie -p 12 --rf -o transcripts.gtf output/SRR5478599.sorted.bam
find . -name *.gtf > list.merged.txt
~/Biotools/gffread/gffread Reference/chr1_pra.gff -T -o Reference/chr1_pra.gtf
~/Biotools/gffcompare/gffcompare -r Reference/chr1_pra.gtf -R -s Reference/chr1_pra.fna -i list.merged.txt
基因水平定量，获得TPM和FPKM的值
stringtie --rf -e -B -A gene_abund.tab -G gffcmp.annotated.gtf -o gene_abund/transcripts.gtf output/SRR5478599.sorted.bam
输出文件 gene_abund.tab有每个基因的表达量值
gene_abund文件夹中的t_data.ctab中有每个转录本的表达量，但是这个表达量只有FPKM的值
这里参考https://www.jianshu.com/p/38c2406367d5
计算count值
这里参考 https://cloud.tencent.com/developer/article/1625215
stringtie的官方文档 https://ccb.jhu.edu/software/stringtie/index.shtml?t=manual
cd gene_abund
python ../prepDE.py -i sample_list.txt -g gene_count_matrix.csv -t transcripts_count_matrix.csv
class code的解释
https://www.jianshu.com/p/5b104830751b
https://www.omicsclass.com/article/234
今天先到这里， 未完待续
网友评论