文献:The Tomato Translational Landscape Revealed by Transcriptome Assembly and Ribosome Profifiling
根据文献,从GEO数据库下载原始测序文件,RNA-seq双端100bp,Ribo-seq单端50bp,两种方式各三个生物学重复。
module load sratoolkit/2.9.6
prefetch --option-file SRR_Acc_List.txt
#下载原始测序数据,ncbi,GSE124962,路径~/ncbi/public/sra/
for i in SRR*;do
fastq-dump --split-3 $i #RNAseq的目录下,RNA双端reads转fq
gzip $i
done
for i in SRR*;do
fastq-dump $i #Riboseq的目录下,ribo单端reads转fq
gzip $i
done
The adaptor sequence AGATCGGAAGAGCACACGTCT was fifirst removed from the Ribo-seq data using FASTX_clipper v0.0.14.文献中去adapter的说明,我用cutadapter。
for i in *fastq.gz;do
out=${i%.fastq.gz}_trim.fastq.gz
"cutadapt -a \"AGATCGGAAGAGCACACGTCT\" -O 5 -o ${out} ${i}"
现在Ribo-seq得到SRR8434774_trim.fastq.gz,SRR8434775_trim.fastq.gz,SRR8434776_trim.fastq.gz,三个文件。
作者用番茄SL2.5的参考基因组和ITAG2.4去掉了测序中的rRNA,tRNA,snRNA等等,用SL3参考基因组ITAG3.2版本注释去出数据中基因组重复序列。将SL2.5和ITAG2.4中对ncRNA的注释下载下来,根据注释提取序列,并且建bowtie2索引。
bedtools getfasta -fi S_lycopersicum_chromosomes.2.50.fa \
-bed ITAG2.4_infernal.gff3 -fo internal.fa
module load Bowtie2/2.4.1
bowtie2-build internal.fa internal
#将ncRNA的序列提取出来,internal.fa,并且建比对索引
将RNA数据比对到ncRNA的序列上,利用samtools将未必对上的提取出来,再用bedtools将bam文件转化成双端fq文件。
for id in *_1.fastq.gz;
do
name=${id%_1.fastq.gz}
one=${name}_1.fastq.gz
two=${name}_2.fastq.gz
out=${name}.sam
bowtie2 -x /public/home/yliang/ly/protocol/tommato-ribo/\
ref/SL2.5/internal -p 10 -1 ./$one -2 ./$two -S ./$out
done
#RNA-seq数据批量比对到small RNA上
module load SAMtools/1.9
for i in *sam;do
out=${i%.sam}.sort.bam
samtools view -@ 10 -bS -f 4 ${i}|\
samtools sort -@ 10 -m 32G -n -o ${out}
done
#将未必对上的read取出来,并且将bam文件排序。
#samtools -f 4指输出没比对上的reads
#samtools sort中默认按照参考gtf文件中的顺序排序,-n加上按照fastq文件的顺序
#ribo和rna的数据这一步都一样
#对于RNA的数据,将bam转化为双端fq,并压缩
module load BEDTools/2.27
for i in *sort.bam;do
fq1=${i%.sort.bam}_sort1.fq
fq2=${i%.sort.bam}_sort2.fq
bamToFastq -i ${i} -fq ${fq1} -fq2 ${fq2}
done
gzip *fq
撒币了,bowtie2有专门输出没比对上的参数,绕了个大圈。对于单端的Riboseq的reads有:
for id in *trim.fastq.gz;
do
name=${id%.fastq.gz}
out=${name}.sam
un=${name}.fq.gz
bowtie2 -x ~/ref/SL2.5/internal -p 10 \
-U ./$id -S ./$out --un-gz ./$un
done
对于双端的RNA-seq数据,有:
module load Bowtie2/2.4.1
for id in *_1.fastq.gz;
do
name=${id%_1.fastq.gz}
one=${name}_1.fastq.gz
two=${name}_2.fastq.gz
out=${name}.sam
un=${name}.fq
bowtie2 -x ~/ref/SL2.5/internal -p 10 \
-1 ./$one -2 ./$two -S ./$out --un-conc ./$un
done
以上步骤去掉了注释的rRNA,tRNA,snRNA,snoRNA。
下一步是去除包含基因组中包含重复序列的reads,下载SL3参考基因组和ITAG3.2对重复序列的注释文件,将重复序列提取出来,建比对索引。
bedtools getfasta -fi S_lycopersicum_chromosomes.3.00.fa \
-bed ITAG3.2_REPET_repeats_agressive.gff -fo repeat.fa
bowtie2-build repeat.fa repeat
与以上相同的方式,通过bowtie2的--un,--un-conc参数把单端Ribo-seq和双端RNA-seq中未必对上的reads给取出来。OK,过滤完毕。